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Prediction of 


The 


Nonlinear 


Pitching and 


Yawing Motion of Svmmetric Missiles’ 


CHARLES H. MURPHY, JR.* 
Ballistic Research Laboratories, Aberdeen Proving Ground 


SUMMARY 


The problem of predicting the combined pitching and vawing 
motion of an aerodynamic configuration which is acted on by 
nonlinear moments has been considered too difficult for a gene ral 
solution. When the problem is restricted to symmetric missiles, 
which may be spinning, relations for the amplitudes and fre 
quencies of the actual motion can be derived by an equivalent 
linearization technique For cubic nonlinearities and relatively 
short trajectories, these relations have been checked by free 
flight tests on the Ballistic Research Laboratories’ spark ranges 

By the use of an “amplitude plane’’ the dependence of the 
character of the motion on initial conditions is simply displayed 
As examples of this technique, amplitude planes are constructed 
for missiles whose moments are cubic functions of the total angle 
of attack 

The ability of the theory to specify possible limit motions is 
described. Theoretical predictions of these limit motions are com 
pared with exact numerical integrations, and the predictions are 


shown to be very accurate 


SYMBOLS} 


aerodynamic coefficients 

pSl/2m) [Cv, — 2Cpn — 

moments of inertia 

amplitude of j7-mede 

dimensionless axial radius of gyration 

dimensionless transverse radius of gyration 

aerodynamic moments about missile-fixed Y, VY, Z 
axes, respectively 

reference length 

pSl/2m) yk.~? Cag 

mass 

I,/I,) (bl/V 


components of the angular velocity of the missile 


gyroscopic spin 


in space along the XY, }, Z axes 


= reference area 


Presented at the Aerodynamics—II Session, Twenty-Fifth 
Annual Meeting, IAS, New York, January 28-31, 1957 

+ This paper is an abridgement of a dissertation submitted to 
the Faculty of Philosophy of the Johns Hopkins University in 
conformity with the requirements for the degree of Doctor of 
Philosophy 

* Deputy Chief, Free Flight Aerodynamics Branch 

t Wherever it is possible these symbols are those defined in 


eference 6 


ov 
{ V/ldt dimensionless distan 
e 0 


yectory 

pSl/2m) |Cvg — ¢ t R Cu, 

components of the linear velocity of the missile 
along the Y, ), Z axes 

magnitude of missile’s velocity 

components of aerodynamic force along tl 
sile-fixed Y, }', Z axes 

angle of attack 

angle of sidesli 

cosine of angle between missile’s axis and trajectory 

sine of angle between missile’s axis and trajectory 

exponential damping of the j-mode 

v+ iw)/V 

air density 

phase of the j-mode 


g2 T ¥ gi ~ ¥ 
perturbation of phase of the j-mode [Eq 
Superscripts 
indicates derivatives with respect to time / 
indicates derivatives with respect to distance 


indicates quantities in nonrotating coordinate system 


> —_ get S pdt 
L.c., & = & 


indicates parameters of basic epicycle without damping 


INTRODUCTION 


N IMPORTANT CURRENT aerodynamic problem is the 
A prediction of missile motion for configurations 
acted on by nonlinear aerodynamic moments. Al 
though linearized theory has proved very successful in 
most cases, there remain important exceptions for which 
the nonlinear terms must be considered. Examples of 
these would include sideward fire of rockets from mov 
ing ships, sideward fire of bullets or missiles from high- 
speed aircraft, and roll-pitch interactions of high-speed 
aircraft. 

Unfortunately, most of the results which have been 
obtained in the field of nonlinear mechanics have been 
restricted to second-order systems and, therefore, are 
not directly applicable to the general sixth-order sys- 


tems which describe missile motion. In this paper we 
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will not consider this general case but will restrict our- 
selves to the intermediate case of a specia] fourth-order 
system. This intermediate case is based on the assump- 
tion of constant drag coefficient and constant spin-to- 
velocity ratio. We will further restrict ourselves to 
missiles with a plane of mirror symmetry and an angle 
of rotational symmetry less than or equal to 120°. 
(These restrictions exclude the important case of a 
high-speed airplane. ) 

Under these conditions special types of nonlinear 
pitching and yawing motion have been studied in refer- 
ences 1-4. By an extension of the Kryloff-Bogoliuboff 
equivalent linearization technique’ the most general 
motion will be considered in this paper. The results of 
this technique can be presented in the simple, graphic 
form of anamplitude plane. The accuracy of the theory 
has been verified by actual free-flight tests and by nu- 
merical integration of the exact equations of motion. 
Some of these comparisons will be presented later in the 


paper. 


EQUATIONS OF MOTION 


Although the linear expansion of the aerodynamic 
force and moment in terms of the dynamic variables 

3, a, B, Pp. q, ", p, 4, * requires the use of 60 stability 
derivatives, the assumed svmmetry reduces this num- 
ber to eleven different nonzero derivatives. Since a 
symmetric missile is frequently given a high spin rate, 
quadratic terms containing p are usually considered. 
These terms are called Magnus terms. For the as- 
sumed symmetry there are only eight different terms 
of this type. This total of 19 distinct coefficients can 
be reduced to seven if only those quantities which 
usually have a measurable effect on the motion are 
retained. 

Under these assumptions the expressions for drag 
and axial spin are’* 


Drag (1/2)pV2SCp (1) 


L = (1/2)pV?SIC,, (pl V) (2) 


As a result of the rotational symmetry the transverse 
components of the aerodynamic force and moment can 
be conveniently expanded in terms of complex vari 


ables: 


—(] 2)pV*SCye (8 + ta) (- 


M+41N = (1/2)p V°S1) [(pl V) Crrpe — 1C uta) X 


, — (¢ + irl 
[8 + ta] + Cw, ' 

Cc i + 1a)/ ipl(B 4 ia) | ; 

l i P OR 
—— V | 
where 
Cx, = —-Cy, = —Cz, 
Cury = Cr a i -_ C, 3 


If the symmetry conditions are relaxed so that a roll moment 
due to deflected control surfaces can be added to Eq. (2), the 


results of the paper will still apply 


Cy oo Ca = ( = oe 

Cy - Cine ( 

Ce one = —C, 3 | 
(The expression following Cy, is the derivative of 8 4 


As has been | 


pointed out by Sacks,’ this derivative is more conven 


ia in a nonrolling coordinate system. 


ient than a derivative in missile-fixed coordinates 
It should be noted that, if the aerodynamic force and 
moment in Eqs. (3) and (4) are nonlinear functions, the 
coefficients are defined by Eqs. (3) and (4) and not by 


differentiation. When these coefficients are functions 


of only a and 8, this definition is unique. This dis 
tinction is important for the nonlinear case. 

When large angles are considered, the exact defin; I 
tion of 8 + 1a becomes awkward, and the use of trans | cori] 
verse components of the velocity seems to be more con- | ¢jaJ] 
venient. We will therefore replace 8 + /a in the above ingt 


equations by a complex quantity € which reduces to it | with 


for small angles and is defined by the equation = 
© | . cond 
= (v+iw)/I oy 
orig! 
According to the above definition, £ is measured ing | of tl 


mnissile-fixed coordinate system. Since nonrolling co- | If th 


ordinates have certain attractive features, we will make | and 
use of a £ which is determined by velocity components | stab! 
. . . = PS {fp dt ° - ~ 
in a nonrolling coordinate system (& = ge’ ”*). Ifthe | mou 
effect of gravity is neglected, the aerodynamic force is 
assumed to be linear, and the mass distribution is 
assumed to be symmetric (J, = /-), the following dif- Th 
ferential equation for € can be derived :* sas 
cn” , , . _— oeth 
Ev + [HT — (7y'/y) — iPl’ — (MW + 7PT)E=0 (6 -_ 
ovmi 
where nl 
Mi 
H = (p.S/ 2m) | ¢ Va 2( a, oe k (¢ Mg T YMG : 
M = (pSl/2m) |yk:-? Cue + (pSl/2m) X \t 
(Cue — Cp) (kh; Ce, + Cpe at. 
SS] 2m)y P : 
a se Aste k ( iat 
1 = ( p.S/ 2m) \¢ Va Up 4 vk, ‘Crsel lamp 
> - ] , \ : 
I = (pl ). CF r Vay 
y = ul (cosine of angle between missile’s axis | jo, 
and trajectory 
yy = Vihar = Vif transverse radius | 
of gy ration 
k, = VI,/mi? = axial radius of gyration the so 
cle | 
and primes denote derivatives with respect to the di 
mensionless arclength s = f (1" /)dt. 
For small angles (y’ = 0, y = 1), constant spin-t 
ws ere 
velocity ratio (P 0), and constant moment coefli 
cients, Eq. (6) can be easily solved: 
g = Ki, e T K c 7 ML, 
ine 11 
where r sta 
issile 
, 
? = pr ¢ # Phe 
- - N 
K = K C enera 
, / = ») > 
ca er - oa ; 
\; + ig 1/2) |-H if $ quit 


V4iM + IP — P? 4 


e oI 6 + 
has bee 
> CONVen 


‘dinates 


force and 


tions, the 
id not by 


functions 


This dis 


ct defini- 
of trans 
nore Con- 
he above 


ices to tt 


ured ina 
ling co 
vill make 
nponents 
). Ifthe 
e force is 
yution 1s 
wing dil- 


ile S axis 


e radius 


spin-to 


it coell 


H 





PREDICTION OF NONLINEAR 
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TION PLANAR MOTION 
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Nonrolling statically stable symmetric missile 


[he pitching and yawing motion 1s, therefore, de 


scribed in a complex plane by the sum of two exponen 
tially damped complex vectors which are rotating with 
ugular velocities ¢,;’. For a statically stable missile 


0) the motion is described 


with zero spin (lJ < 0, P 
by a damped ellipse (see Fig. 1). For special initial 


Ky K3), 


origin, the ellipse degenerates into a line, and the axis 


conditions the motion goes through the 
of the missile is confined to a plane (planar motion). 
If the missile is allowed to spin, the ellipse will precess 
and become an epicycle (see Fig. 2). A statically un- 
stable missile which is spin stabilized has an epicyclic 


motion with the nodes on the inside. 


SOLUTION OF NONLINEAR EQUATION 


The most common nonlinearity associated with a 
symmetric missile is a dependence of the aerodynamic 
coeflicients on the sine of the total angle of attack. 
Symmetry requires that they be even functions of this 
It 


namic moment coefficients in Eq. (6) 


uantity. is, therefore, assumed that the aerody 
are functions of 


c 


At this point we impose the reasonable condition that 
for “small” nonlinearities the solution of the nonlinear 
equations can be approximated by an epicycle without 
lamping. The approximate solution is then improved 
by a perturbation method. If Eq. (6) is written in the 


lorm 


Ps’ (0) € —[H (n?) — (y' y) Jf’ 4 
AI (n — J/(0 + iPT(n |; (S 


the solution of the left side set equal to zero is an epi- 


cycle without damping. 


&€&= Kye" + Koo 9) 
where 6; = do + 6,’ s 
6) = (1/2) [P+ VP? — 4M(0 
r 117(0) > 0 


he inequality, P? — 4.1/(0) > 0, is always satisfied 


lor statically stable missiles. For a statically unstable 
missile it is the criterion for proper spin stabilization. 

[he perturbation procedure which will be used is a 
generalization of the Kryloff Bogoliuboff method and 
is quite similar to the method of variation of parameters 
lor the solution of inhomogeneous linear equations. 


P | 


scout 


NG AND YAWING MOTION $7: 


We assume a solution of the form 


£ K ¢ - Ko ’ 10 
where Aj(s) and y,(s) are unknown functions of 
Differentiating Eq. (10), 

g’ 1d,’ Aye TOY be ide K v6 A 
1, Ay t — A T y ‘A 1] 


The sum of the last two terms is set equal to zero, and 


£’ is differentiated. 


(A,’ T Wy’ Kye —? + 
(Ko! + tWo’Ko) ¢ . 0 12 
&” —o"*Kye re) _ ¢."*Kx . o (K 
iW,’ Kye’ si T oo (K . l fj . 13 
If Eqs. (10), (11), and (13) are substituted in Eg. (S), 
then 
id’ (Ky! + 1Wy’Ky)e ht) 4 56."(Ko’ 4 
tWe'Ko)e'" Tr" —|I1(n y yié 
[.\/(m W(O) + iPT (n°) |g 14) 
where 
n é = A,? + A.* + 2A,A> cos 
Y ¥i<- n 
QO = Qo T Yr —_— es — V1 
and ¢ is given by Eq. (10). Ae’ + iWe’A» can be elimi 
nated between Eqs. (12) and (14) to yield the following 


equation for the functions A, and ¥ 


(Ay) Ay) + ty’ Ky, Ko, @ 5 
where 
Ay, AK , @ QO) oO T/ n 
(y’/y)] [oi’ + de’ (Ke Ay) ,;PT(n 
(nm (0) ]; [1 KA» Kye }) 
Eq. (15) is very difficult to solve because of the ex 
plicit presence of the independent variable in ¢. Our 


basic assumption that the motion can be approximated 
by an epicycle without damping implies that the de 
f the 
This equation 


are such that the character « 


rivatives in Eq. (15 
motion changes little over a cycle of @ 
must describe changes in phase and amplitude over a 
We therefore a 


a cycle of @ and separate into real and imaginary parts. 


number of cycles. verage Eq. (15) over 





+ ‘kK 
K 
~~ ‘ 
7 
r “2 

a at a A eee ee a eee 3 
| 
4 
FIG. 20 FIG. 2b 
Fic. 2. Rolling statically stable symmetric missile 
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Ki'/K, = (A, Ke) 
ag te | 
= ? ) LT (n?) [di’ + 
G2) J0 


. , 
23 (hi — 


o»'(K» Ky) cos ¢] — PT(n? x| 
[1 + (K»/K,) cos |i d@ 


wy 


—] 
Wy! = : 11M (m?) —(0)] X 
YQ 0 


29 (d,’ — Po je 
[1 + (K2/K,) cos d] + (y’ vy) X 
go'(K2/K,) sin df dd (17) 


Similar equations for A, and y. can be obtained by 
interchanging the subscripts | and 2 in the above equa- 


tions. 


APPLICATION OF THEORY TO SHORT TRAJECTORIES 


If we consider short trajectories over which the modal 
amplitudes (A,’s) do not change very much, the solu- 
tion obtained in the preceding section can be easily 
described. The motion is a damped epicycle with 
frequencies ¢;/ = ¢,’ + yw,’ and damping exponents 
\,;. Both the frequencies and the damping exponents 
are, however, functions of A; and A». 

An interesting application of the theory has been 
made in the data analysis of test firings on the Ballistic 
Research Laboratories’ aerodynamics ranges. For 
cubic static moment and small angles, 


M = My + Mon? (18) 
y = 1 — (1/2)r (19) 


Substituting these equations in the equations for the 
frequencies, we have the following equations: 


oy’ = Gy’ — [12(Ki? + 2K2”)/ (gi! — d2’)] — 


(go’Ke?/2) (20 
p 7 ane , t 2U) 
oe’ = do’ oe | VWo( Ko? + SK) (do’ — o;') | _— 
(¢,'K\ 2) 


where my = G/2) |P + Vv P? -- LM 


If these equations are multiplied and the A,;* terms are 


neglected, 

gr!’ -de' 1 + [(Ki2 + Ke?) /2]} = My + Mon2 (21) 
where 

ne? = Ky? + K.? + [(gi'Ki? — 2'K2")/ (br — o2’)] 


From a group of flight tests of the same configuration 
at the same Mach Number, values of A; and ¢,’ can be 
obtained. These values can be fitted by Eq. (21) to 
yield 1/) and \/,. This has been done for a number of 
programs.* The results possessed excellent internal 
consistency and compared well with direct wind-tunnel 
measurements. Similar relations have been used to 
measure cubic Magnus moments. 


THE AMPLITUDE PLANE 


Although the case of short trajectories is of some im- 
portance, the primary interest of a designer is in the 
behavior of the missile over long trajectories. As 
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long as the motion predicted by the linearized equatjy a 
is gyroscopically stable [P*? — 4.17 (0) > OJ], it can 
seen from Eqs. (16) and (17) that the only effect of .| 
mildly nonlinear static moment .\/ is to change the fr 
quency of the motion. Unfortunately, strong pn 
linearities in ./ would invalidate the perturbat; 
process and, hence, are not considered in this paper 
Strong nonlinearities in the other moments affect {J 
damping and can be considered if the predicted dany 
ing over a cycle 1s reasonably small. It is this dampin; 
or undamping of the motion over a large number 
cycles which we will consider in detail. 
In order to do this we will study the following firs 
order nonlinear equations: 
| 
(K,*)'/ Ki? = 2\(Ki?, Ke ne ¢ 


(Ko")'/K 2r.(K,", Ko? 2 


where the \; (A,’, Ay”) are given by Eq. (16) and th 
corresponding equation for the 2-mode. These equ 
tions are now divided to eliminate the independer 


variable. 


d (K»?) d (A?) = K o“\ ts, K-’ Kk, “\y Kk, 7” Z 


Eq. (2+) has been extensively studied by Poincar 
Its essential behavior is determined by its propertie 
in the vicinity of points for which the numerator a1 
denominator both vanish (singular points). 

Solution curves of Eq. (24) are trajectories in th 
K,’-K.” plane. This plane, which will be called th 
amplitude plane, is of fundamental importance fo] 
stability studies of missiles acted on by nonlinear? | 
moments. ; 

For the linear case, the \,'s are constants, and tl 
origin is the only singular point. If the spin is zer 
A; = de, and all the trajectories go either toward ot 
away from the origin (see Fig. 3a). The particular 
curve which describes a given motion is specified by th 
initial values of A; and direction of motion along t! 
curve by the sign of \;. For positive \,’s the motic 
undamps away from the origin. The damping @ 
ponents of spinning missiles are not necessarily equ: 
and may differ in sign. The origin can, therefore | 
be either a nodal point (Fig. 4a) or a saddle point (Fig. J 
fb). All of these figures display an important chat 
acteristic of linear systems—i.e., the type of motion 
independent of initial conditions. 

Amplitude planes will now be constructed for tw} 3 
simple nonlinearities. In the first example, the effect] 
of a nonlinear damping on the motion of a nonspinnins 
missile will be considered while, in the second, the effec! 
of a nonlinear Magnus moment on the motion of 
spinning missile is described. In both cases the angle] 
will be assumed small so that the geometric nonlineat* | 
ties can be neglected (y’ = 0, y = 1). 

Example 1: 

=. ay + a2 7 


Cuy + Cua 


For zero spin, ¢:’ = —d¢»’; and the equations for tht] 
damping exponent \, take on the form 
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(KY, Ke? (1/40) fo 0) y, In Fig...5, the origin is a node, and the dep nden 
Jo of the type of motion on initial conditions is once mor 
[1 + (A2/ Ki) cos @] dp (25) shown. Figs. 6 and 7, which illustrate amplitud | 
- planes for which the origin is a saddle, show limit m 
Ne (K;?, Ke?) = (1 +427) | II(n?) X tions. In Fig. 6a, the nodal singularity on the A 
ibs ee ; predicts a limit circular motion of the missile’s ayjs | 
[1 + (Ki K:) cos g] dp (26) about the trajectory. In Fig. 6b, the spiral singularjt 
where n = Ki? + K.? + 2K,K2 cos o predicts a limit epicycle. Fig. 7 shows the transitig, 
If the damping derivatives are linear functions of 7”, ae ni snd ones mami ste ; i . 
* d : i. é g g € plane 
E(g") as the stmpte form In summary, the important problem of stability fo 
H(n?) = Hy + He 9? (27) nonlinear motion can be reduced to the construction 
the appropriate amplitude plane. The initial con 
where tions for which the motion will remain bounded ¢ay 
Hy = (pSI/2m) [Cveq — 2Cpd — kt~? ao] then be easily identified. 
HT = —(pSl/2m) k,~? ar 
NUMERICAL CHECK OF THE THEORY 
\(K17, K.?) = —(1/2) [Ho + 2K,’] (28) a 
Ihe examples in the preceding section show that A 
ho (K,2, Ke?) = —(1/2) [Ay + TK?) (29) rapid calculation of the dependence of the motion oy | mot 
The amplitude plane for the interesting case for which mens rege oe 8 oe role — = 
the damping changes sign (//)//, < 0) is shown in Fig. —_ oie — —— oo oe > this 
3b. The lines of zero damping are given by Eqs. (28) mae ne paeieneraneaien hay meas oo val 
and (29) set equal to zero and are shown in the Figure mane onamatanse may compe — meen We ae oF 
together with the four singular points at [0, OJ, a may elagiaeogianee ig — “A the 
(0, —(Ho/H»)], [—(y/ He), 0], and [—(U/Ib), in actual missile flights. It remains, however, to d met 
— (Hy/ Hs) |. termine the quantitative accuracy of the theory. 
If the origin is first assumed to be a stable node, — _ — pore ms ae " — abet : 
this Figure shows the characteristic nonlinear depend- en ~ 3 ~ om asain PPOveng Ground 
ence of the type of motion on initial conditions. For wae Se: Ene, henge ee ae SE ee om 
initial amplitudes which are inside the square formed nig ——. seni — a reper 1 sy 
by the axes and the zero damping lines, the motion will rinlie paggtr an ree tt hg ep gs leer a M : 
damp to the origin; for amplitudes outside, the motion pane wou ao ‘aguaeae nid mai une. - 
will grow. amplitude plane for this case 1s shown in Fig. 8. | Rock 
Under the reverse assumption that the origin is an ro, a meres ora a 
unstable node, the other node becomes a stable one, — we — . decag nett oe ~ y 
and all trajectories lead to this limiting value. This esti ree jee ne Se Se a vaemes ot K 
existence of a limit motion is another characteristic of 7 a i ge pre rien eres es ro 
nonlinear systems. It is interesting to note that, if we pants a ¥, ane Re nage = nen ae eee “ 
had restricted ourselves to planar motion, Eq. (6) creased were denoted by U (unstable), those nh 
would have reduced to the well-known van der Pol — —_ eam «ie ere -~ _— 4 — 
equation, and the singularity at [—(/J) 7/2), —(U1) + — - or —_ es Se 
Hi.) 1 ia the Temit cycle predicted for that equation. since the modal amplitudes do not completely deter 
Example 2: mine the initial conditions (&’, &)), the initial cond 
tions were selected for both an initial maximum and at 
P, H constant; Cy. = bo + b. 9 initial minimum of £. The theory makes no distinction 
; ; between these cases. As can be seen from Table 1, the 
The equations for the damping exponents reduce to ee ee : ‘ 
predictions of the theory are quite good. 
Ar = Aw + Ave (Ay? + 2A2’) Next, the existence of the predicted limit motion was 
Ne = Dew + Aes (Ke? + 2K,") investigated. This was done by selecting initial con 
: ; alias ditions corresponding to amplitudes which lie to the = 
where left of the separatrix and integrating for 30 periods 0! 
oe Oe ~~ Nn’. The motion near the end of each run was analyze 
eo = (H de’ — PT.)/(be’ — dh’) and found in all cases to be epicyclic with zeTO damping. 
i se colle we wh = 6 rhe values of the amplitudes checked to four cecum 
Ty = (pSI/2m) (Cue — Co + Ra? bo] places with the predicted values of A, ad 0.0268, kK; : 
T, = (pSl/2m) ka? be 0.0784. I'wo of these runs are shown in Fig. S, " 
? the time history of Run 1 is shown in Fig. 9. Tht 
The quadratic Magnus moment coefficient which cor- curves in Fig. 9 were obtained from Eqs. (22) and (2 


responds to a cubic Magnus moment can give rise to while the circled points came from integration of Eq 


several interesting amplitude planes (see Figs. 5-7). (6). 
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TABLE | R 
Initial Initial 2-2 2 
107A Maximum Minimum . : - P 
HO 0. SOO S Ss 106 | 
600 O.825 >» 5 
HOO QO 850 S cS K \ 
600 0.875 S S » Ky, 
» 600 0.900 N | 1044 
600 0.925 U l 
Predicted 10°A,? O.895 p f P 
200) 0) 150 S S ooh - 
1 200 0.475 S S _ 
1.200 0.500 S S » 
1 200 ().525 S S \s ah eatrcetaietl 
200) 0.559 S I > an ene corre } 
P00) 0.575 [ | ae = 
D1) 0.600 I I eriods of 7 
Predicted 107K 0.555 Fic. 9. K,2 and K>2? vs. periods of 
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SUMMARY needed for a fundamental study. It has been suggested 
Studies of transition over flat plate at Mach Number 1.76 that a brief account of these results would be useful t 0 
were carried out using a hot-wire anemometer as one of the other investigators in the light of our present under 
principal tool Phe nature and measurements of free-stream standing of the problem, at least in a qualitative man- ' 
disturbances a Ipersouic speed ire inalyzed rhe experi ner. In particular, it should throw some light on th: 
mettal results are interpreted in the light of present overall in > a : me ees 
: . possible influence of free-stream ‘‘turbulence’’and other 
formation on transition at supersonic speeds and conclusions as ; e I 
to further fruitful experime its are drawn tunnel disturbances on the differences in transitior corn 
Reynolds Numbers observed in different tunnels. For Mac 
SYMBOLS flat plates in streams near Mach Number 2, reported valu 
ree : transition Reynolds Numbers range from approximately line : 
db decibel level of souid intensity with the pressure of he 5 mes ; _— late Mne 4 
& xX 1074 dvnes/em:? as reference level 2 log, one million to more than three millions. smal 
Ap’) + 74 chan 
V/ Mach Number EXPERIMENTAL CONDITIONS face 
p pressure, dynes /em = chan 
Re Revnolds Number, based oa free-stream conditions and The test section of the tunnel expands slowly to con 
: ; stagt 
length indicated by subscript pensate for the growth of the wall boundary laver. It is , - 
eat solute te er — or: . 52 1: . the t 
f mean absolute temperature 7.27 in. wide and 10.93 in. high at the leading edge of th 
verturbation in axial velocity compo.teut ; : _ : ber 
,' ; flat plate, where the thickness 6 and the displacement 
l mean axial velocity not s 
oe ic Ae at thie diaieics Pt ae en , 7 
vertical velocity component thickness 5* of the turbulent sid all boundary lave saa 
\ ixial distance from leading edge of plate, in were approximately 0.6 and 0.12 in. In order to mini d 
: ; - ° : ind | 
vertical distance from plate, in mize the interaction with the sidewall boundary layer 
: IC 
1 dis e from ceaterline of plat : : , 0.0] 
latevat distance from cestertn per the plate extended only about halfway into the layer at 
6 = boundary-layer thickness a erie layer 
; 5 ; ; : the leading edge, where it was 6.75 in ide, and wv , 
6 = displacement thickness of boundary layer 
: 3 . . ar > ) } re 15 lownstrean 1) t¢ N 
A perturbation of the qua itity which follows attached to the wall 15 in. downstream. The pl 
A momentum thickness of boundary lave 0.25 in. thick with a leading-edge angle of 11.2 set 
Sisbscriprs Fig. | The higher pressure from the slanted lower sur It 
0 free-stream value ia settling chamber face ‘leaked’ through the sidewall boundary layer over exce 
| = free-stream superso vic valuc a distance of 1.26 in. and was followed by an expansiot dom 


Superscripts 


wave of approximately the same strength and lengtl 





boun 


root mean square value of fluctuation The adverse pressure gradient on the top surface was | in ch 
nohec r ‘ : oe 1 
a thus kept within two oblique strips (g¢ ~ 35°) from the | exam 
as 
if 
INTRODUCTION tive 1 
~ -<«) e » 7 S I . 
—_— 1952 and 1955 a series of experiments on o. ~ 
rine ° ° , v b 
transition on a flat plate were performed in ‘i he P, oun 
- — ae ; Si rate as 4 - a travel 
Johns Hopkins University continuous supersonic wind 7 saa ciligs calainia jos tes . malt 
tunnel at the free-stream Mach Number of 1.76. The ee y “ntry 
ne ; ; . ‘ i / nounc 
distinctive feature of the program was the usage of a 3| SG unc 
: amok a * the tr: 
supersonic hot-wire anemometer. It was concluded at -z= 10 — 
‘ . F ‘ - P =7 2-15 probe 
the time that, with the available equipment, the ind1 2| seine 
« - Pre 4 [ Wire O 
vidual measurement and control of the many causes of ’ 
“43 , very 
transition and the measurement of the corrcsponding aes . 7% 
, 8) —») > 
: : ° ° i 5 nal oO 
disturbances in the boundary layer fell short of that een ee edge of ” 
ntermiftent bursts 
Received October 22, 1956 es C = rar we? 1 t Si 
+ This research was supported jointly by the Office of Nava! 2 aren aera ey er A gr 
‘ ‘ . a ae : +) 2 4 6 Rg , 
Research, U.S. Navy, and by the Air Office of Scientific Research, x vibrat 
U.S. Air Fore The author is greatly indebted for many-sided Fic. 1 Edge of boundary layer and wall pitot t1 A] 
help and cooperation on the part of his colleagues and assistants, transition region = pitot pressure with probe pressed lave 
in particular Dr. F. Clauser, Dr. L. S$. G. Kovasznay, Mr. L. T plate 2 alias free-stream total pressure distance — ae vibrat 
4 ts ‘ line; 1 beginning of mild adverse pressure gradient 1! 
Miller, Mrs. F. Baake Chenoweth, and Mrs. P. Clarken Mellroy corner of plate; | = beginning of compensating favorable pres wav 
Research Scientist, Department of Aeronautics sure gradient; Reynolds Number per in. ~.387 X 10° 
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Mean velocity profiles ahead and near end of transition 
region 


corners and corresponded to a maximum (gradual) 
Mach Number change from the unaffected forward 
value of 1.72 at x = 5.0 to 1.69 at x = 6.25 on the center- 
linez = 0. For most of the tests the plate was kept at a 
small negative angle of attack in order to decrease the 
chance of a disturbing separated bubble on the top sur- 
face.*> Since there was some variation of the settling 
chamber turbulence level over the available range of 
stagnation pressures (+10 per cent of atmospheric), 
the tests were run at an essentially fixed Reynolds Num- 
namely, 0.37 & 10°. Continuous (but 
vertical and horizontal traverses with 


ber per inch 
not simultaneous 
standard JHU (Johns Hopkins University) hot wires! 
and conventional circular Pitot tubes (0.016 in. OD, 
0.010 1D) provided the observed variations of boundary 


layer characteristics. 


TRANSITION PATTERNS 


It does not appear to be generally appreciated that 
except in the regions where Tollmien-Schlichting waves 
dominate) the traces from a hot wire inside a laminar 
boundary layer differ primarily in amplitude but not 
in character from those in a turbulent layer (see, for 
example, Fig. 13 of reference 2 and Fig. 11 of reference 
17). A clear distinction and, therefore, a simple posi- 
tive identification of the transition region can, however, 
be made from the character of the entry into the 
boundary layer from the free stream during vertical 
traversing with the hot wire, especially at high speeds. 
Entry into the laminar layer was characterized by pro- 
nounced low-frequency (90-120 c.p.s.) undulations of 
the trace corresponding to small natural vibrations of the 
probe-supporting mechanism; vertical motion of the 
Wire on the order of 0.001 in. across the ‘“‘edge’’ of the 
very thin layer corresponded to effective velocity sig- 
The outer edge of the 


nals of 3 ft. see. and higher.+ 


TSimilar generation of signal by lateral motion across a shear 
gradient was employed to determine the natural frequency of 
vibration of the flat plate 
A low-speed jet in the otherwise quiescent tunnel created a shear 


another potential transition promote 


iver over the plate When the plate was struck sharply, its 
Vibrations moved the layer across a hot wire, giving rise to regu 
‘T Waves between 180 and 220 eps, of no consequence to super 


transition 


turbulent layer is detected from the unmistakable high 
intensity short-duration spikes corresponding to the 
difference of fluctuations in the free stream and in the 
turbulent crests of the layer (see Fig. 16 of reference 1 
Sample outlines of the edges of the layers so obtained 
are seen in Fig. 1. 

Fig. 1 that the 
originate well within the laminar layer but that they 


suggests initial turbulent bursts 
grow very rapidly outward, at lateral speeds of 30 to 60 
ft. sec. Limitations on spatial and frequency resolution 
due to the high speeds! prevented more detailed study 
of the structure and of the statistical distribution of 
these spots, which could well correspond to those of 
Emmons.* * From continuity considerations, the up- 
stream edge of the turbulent fluid should radiate a local 
narrow rarefaction wave. However, the occurrence of 
the high lateral speeds of the turbulent crests above the 
extrapolated edge of the laminar layer suggests that 
actual motion of mass is involved and not merely a 
wave propagation of velocity in an unstable layer 


Hence sound waves should be 


generated whether the relative horizontal speed of the 


strong compressive 
turbulent lumps with respect to the free stream is sub 
sonic or supersonic. Horizontal traverses with the hot 
wire well above the boundary layer indeed picked up the 
radiated sound as a substantial but gradual rise in a.c. 
signal downstream of the Mach line originating near 
x = 6.5, Analysis’ of the incremental signal indicated 
that the wave fronts of the sound pulses do not exhibit 
any dominant direction, even though as a whole they 


move outward along free-stream Mach lines. Intensi 


ties on the order of 140 db. or more (based on 2 X 10~4 
dynes cm.” as zero pressure level), that is, r.m.s. pres- 
sures of 0.9 per cent of local free-stream pressure were 
inferred. The fine structure of the field (which at higher 
speeds includes discontinuities observable on shadow- 
grams)®** should contribute substantially to the total 
sound energy but was obscured and partially suppressed 


by the frequency limitations. 


COMPARISON OF TRANSITION REYNOLDS NUMBERS 


Repeated hot-wire breakage close to the plate (es- 
pecially in the violent region between x 6.5 and 7.5 
and local flow separation’ when the probe approached 
the wall farther upstream prevented exploration of the 
edge of bursts ahead of the x = 6.S station. Reasonable 
extrapolation of the hot-wire measurements from Fig. 1 
would lead to the following Reynolds Numbers for the 


startt of the transition region: 

t Since one could not reliably distinguish between signals from 
small laminar lacunae imbedded in turbulent regions and signals 
from the variation in turbulent velocities themselves at these high 
speeds the hot wire gave a poor indication of the end of the 
In particular, the criterion of maximum wall 
x 10° while the 


transition region 
pitot pressure indicated an ending at Re, ~ 3.2 
a.c. hot-wire traverse of the layer at the point showed little 
resemblance to the structure of the developed turbulent layers at 
the wall (reference 1, Figs. 16 and 17 


profile could not be fitted into the customary power-law shapes 


, and the mean velocity 


) 


see Fig. 2 
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Re; ~ 1.6 X 104, 
Re, ~ 1.7 X 103 


Ke, ~~ Z.a X 10, 
Res ~ 7.6 X& 103, 


Comparison should be made with the highest reported 
transition Reynolds Number on flat surfaces at the 
same Reynolds Number per inch and comparably “‘low”’ 
supersonic Mach Numbers. Coles’ data’ (see Fig. 18, 
Re/in. ~ .37 X 10% at M = 1.97) give Re, ~ 3.3 X 10 
when judged by the rise in wall static pressure and 
Re, ~ 2.6 X 10° when judged by the minimum of shear- 
ing stress,** while Shoulberg et al. (Fig. + of reference 
13,'T Rein. ~ 0.3 to 0.46 X 10°, J ~ 1.98) indicate 
Re, ~ 3.3 X 10° based on the minimum of the tem- 
perature recovery factor. The minimum shearing stress 
criterion should correspond to that of minimum wall 
pitot pressure (see Fig. 1) which on the centerline gave 
Re, ~ 2.1 X 10°. In view of the known trends®: !° of 
transition Ke with Mach Number, the JHU transition 
would appear to occur somewhat early. 

The fundamental fact which obscures our under- 
standing of transition lies in the existence of large num- 
ber of its causes! (promoters) whose effect cumu- 
lates and varies with local boundary-layer velocity and 
temperature distribution, Reynolds Number, Mach 
Number, etc. Often enough one or two causes domi- 
nate and mask the effect of the others until they are 
sufficiently reduced in strength so that other dominant 
causes may emerge. It is therefore not surprising that 
different testing facilities indicate contradictory trends 
as a given parameter is varied. For experimental clari- 
fication it seems essential to be able both to charac 
terize or measure the causes and to control the corre- 
individually 


nondimensional parameters 


Furthermore, studies$ of the cu- 


sponding 
within wide limits. 
mulating response (prior to transition) of the laminar 
layer to the controlled promoters should furnish clues 
The back- 
ground of even partial theories!’ '* is of course an 
invaluable help for the design of the experiment and 


as to the mechanism(s) of transition itself. 


the correlation of the data. 


THE EFFECT OF VARIATIONS OF DISTURBANCES 


The specific difficulties with supersonic experiments 
center around the problems of the detailed control of 
model and wind-tunnel conditions and of the spatial 
and timewise resolution of the measuring instruments 
(which include probe interference, vibrations, noise, 


etc.) 


** Need for caution when interpreting experiments only in 
terms of location of transition is underlined by the fact that, in 
the same tunnel, with essentially the same plate as Coles worked 
in 1952, reference 25, p. 7, reports in 1956 the beginning of transi 
tion as sensed by a hot wire at Re, ~ 0.9 X 10°, whereas Coles’ 
data interpolated to the corresponding Mach Number and 
Re/in. would indicate Re, of 1.7 to 2.1 X 10° 

'T Reference 13 exhibits a trend of transition Re with Re/in. at 
variance with that generally reported” © for models without 
pressure gradients, especially with reference 9 

§ Such as attempted at low speeds by Schubauer and Skram 
stad,'? Schubauer and Klebanoff,? Bennett,? Liepmann®*® and 
others and at high speeds by Laufer and Vrebalovich * 
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Gross changes in those promoters which were likely 


to affect the laminar layer and the transition pattery St 
as described were tried in order to bracket the effects th 
The distributed roughness of the plate increased grad w 
ually through the experiments and was diminished st 
abruptly by careful repolishing without any palpabk bu 
effects. Trippers near the leading edge, especial = 
those that provided a length of separated boundary | . 
layer (such as two bands of tape 0.014 in. thick witha | jo) 
gap of 0.35 in.) did move the transition pattern for vn 
ward—.e., it became dominant at this relatively loy by 
Mach Number.j The actual location of the transitioy vel 
occurred 0.5 to 3 in. downstreamf of the reattachment the 
but seemed to depend rather strongly on free-strean by 
disturbances (cf. page 6 of reference 21). The uni of s 
formity and sharpness (approximately 0.001 in. in diam ( 
eter) of the leading edge deteriorated with time and | seg; 
was restored without appreciable shift in transition the 
This lack of sensitivity does not contradict our present | for 
knowledge’ of a range of influence of leading-edg con 
thickness because of the relatively small Mach Number | 9.2 
and small thickness changes in the experiment. Span | od 
wise nonuniformities of the leading edge have bee sect 
known to promote earlier transition on wedges pre 
sumably by generating streamwise vortices but aj bule 
parently did not constitute strong enough disturbances to cl 
in the present case. | the. 
Furthermore, no measurable shift of transition r ing] 
sulted when the axial component of the settling chamber | et al 


turbulence+{} was raised from 0.74 per cent to 4.6 per | dime 
cent of the mean velocity (10.2 ft. sec.) by installing ntl 
special grid (0.75 in. square bars, + in. on center pF tl 

' 


in. upstream of the beginning of the nozzle. Tempera 
ture fluctuations in the settling chamber were on the | low 


order of the noise of the equipment—negligible. Th .cros 


settling chamber turbulence was also raised by placing | old 
a horizontal cylindrical bar, | in. in diameter, 12 1 | 
upstream of the beginning of the nozzle and in line with | \ 
the leading edge of the plate. The turbulent patter strea 
downstream of such rods, including the probabilit in 1%) 
distribution and spectrum of «’, was known from | press 
Roshko’s studies.*” The turbulent fields observable 1 speed 
the settling chamber are, of course, strongly modifi¢ Iree-s 
during the passage through the nozzle. These changes | Some 
can be estimated to some extent by relying on linea! re} 
theories of Ribner and Tucker” and on their comparis prob 
with Uberoi’s experiments.** Some limited measut lence 
ments of turbulence in the supersonic test section ar 
also available as will be discussed in the section on Ire FREI 
stream disturbances in supersonic tunnels in general \t 
1eW¢ 

+ Trippers become less effective as Mach Number inet aa 
and it is probably more than coincidence that there 1s soul 
sponding very large increase in the stability of free bound temp 
lavers as seen in Fig. 8 of reference 19 and predicted in refet Irrot 
35 

t As in reference 18 (.1/; = 5.8) lateral contaminati . 
were narrower than at low speeds,* namely 5° to 7°, increasili se 
slowly downstream —— 

ti These and the simultaneous supersonic measurem ATLOW 
performed by P. Clarken and L.S. G. Kovasznay entre 


ere likely 
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The circular bar (Reynolds Number ~ _ 5,400, 
Strouhal Number*" ~ 0.21) was shoved in and out of 
the settling chamber as the heated wire monitored the 


supersonic laminar layer on the plate. The comparison 


established beyond doubt that, with the rod in, tur 


bulent spots, ‘‘lasting’’ less than a millisecond, originated 


UU 


somewhere upstream of x = 3” with an average fre 
quency of approximately one in twenty seconds—+.e., 
not often enough to influence the mean location of 


transition. Presumably, these spots were brought about 


by the high tail end of the distribution function of the 


locit not unlike 


velocity 


the spots of Schubauer and Klebanoff 


fluctuations (reference 20, Fig. 19 
were generated 


by high velocity disturbances associated with discharge 


| 


7 ark 
ol sparks 


Otherwise, the lack of sensitivity of the transition 


region on the plate to changes in stream turbulence in 
he settling chamber contrasts with the effects reported 


t 
t 
for 10° cones by Evvard et al.'* at JJ, = 3.12, and for 5 
cones by Laufer and Marte!’ at Mach Numbers below 


22. However, Laufer and Marte found no effect of 


settling chamber turbulence on transition in the test 


section at .W/; = 2.55 and higher. 


slender cones are more sensitive to free-stream tur 


It is possible that 


bulence with high energy slow eddies being equivalent 
to changes in angle of attack [this might be inferred from 
the observed rapid forward shift of transition when the 
see Witt,’ Fig. 9, and Jedlicka, 


presumably because of lower stability of three 


ingle ot attack changes 
et al 
dimensional boundary-layer flows**]; (2) the transition 
JHI 


lree-stream 


plate was dominated by some other cause 
the 


mum disturbance level in the supersonic section was 


on the 
than 


fluctuations; (3 actual maxi 


lower in the JHU tunnel because of the high speed ratio 
cross the 150 
old observed in reference 10 at higher Mach Numbers, 


i a Q-- 
> +) 


nozzle and remained below the thresh 


\ number of measurements were made of the free- 
stream disturbances at the time of the experiments and 


in 1%)9 but these must be interpreted in terms of our 


present understanding of hot-wire response at high 
speeds lo appreciate the situation let us consider 
iree-stream disturbances in supersonic wind tunnels in 


some detail, especially since many other discrepancies 


1 reported measurements have been ‘‘explained”’ as 
robably due to differences in free-stream ‘‘turbu 
lence See also reference 3°. 
FREE-STREAM DISTURBANCES IN SUPERSONIC TUNNELS 
\t any given point fluctuations in a gas can be 
viewed! ° in the first approximation as a superposition 
three modes of covarying physical properties 
sound’ mode (variation of pressure, density, and 
temperature at constant entropy and oi the coupled 
irrotational velocity field); (2) ‘‘entropy”’ mode (varia- 
tion of entropy, density, and temperature at constant 
pressure ind (3 “vorticity” mode (variation of the 


solenoidal component of the velocity field which is 
snown as “turbulence” at incompressible speeds). The 
entropy and vorticity modes are essentially convected 
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along streamlines (see Fig. 3) so that in a supersonic 


tunnel they must be traceable (with distortion) to 
conditions in the settling chamber (or be generated dur 
ing the acceleration through the nozzle, a second-order 
effect according to present beliefs The sound disturb 
ances can travel across streamlines so that they may 
come from the settling chamber? and from the bound 
the test that 


are high enough in any of the three modes may promote 


aries of section. Disturbance levels 


transition, although the actual coupling mechanisms 
are not known. 

For sound fluctuations originating at 
the wall 


convenience, 
see Fig. 3) can be classified into four types 
(a) radiation from nascent turbulence as described 
radiation from developed turbulent bound 
layers 24 c) diffraction 


and scattering of otherwise steady pressure gradients 


earlier, (b 


references 2: 


ary lsee 
and shock waves (as generated by nozzle contours, un 
intended waviness or roughness, models, supports, etc 

through the turbulent boundary layer, and (d) radiation 
from unsteady wall vibrations caused by pressure fluctu 
ations in the boundary layer or by the loads on the dif 


fusor, associated with the unsteadiness of the terminal 


shock wave.** While the hot-wire anemometer is 

+ In low-speed tunnels sounds from the various parts of the 
system are critical for low-turbulence operations 2 In super 
sonie tunnels the noise level has to exceed 120 db or more to be of 
consequence relative to the high free-stream speeds—1.e., 1,500 
ft. sec. in the present case rhe amplification of the sound due 
to the contraction could not exceed 6 db Sound levels on the 
order of 114 db were definitely not reached in the settling chamber 
as the author can testify from personal experience. As the wind 


tunnel was going up in speed, the noise in the settling chamber did 


rise to very uncomfortable heights until the speed of sound was 
reached in the nozzle and cut off the sound sources in the diffusor 
rhereafter, a pleasant though breezy atmosphere prevailed wit! 
the noise level estimated at no more than 70-SO db —of no conse 
quence for disturbances in the supersonic test sect 
There may be substantial generation of sound when high 
intensity entropy spots are swept through the nozzle, as in Jet 
engines and rockets 
\ssessment of the magnitude of sound of category (d) along 


the lines of Ribner® and Corcos-Liepmann™ indicates that It 1s 





unlikely to reach the intensity of the sound of category (« i 
120-130 decibels \ similar conclusi follows from the com 
parison of the wall sound radiated into the tunnel with its 
counterpart, sensed in the room around the tunnel. Since the 


» the skin of aircraft) are designed to 


tunnel walls (in contrast t 

minimal deflections, this conclusion is not surprising. Other 
factors which limit the inward wall radiation are (1) the very 
low air densitv near the wall, (2) the relatively high density of 
the wall material, and (3) the relative thinness of the boundary 


laver which governs the scale of the pressure fluctuations 
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sensitive to all of the fluctuation modes and cate- 
gories of sound, the combined electrical signal can- 
not be satisfactorily decomposed into the separate 
contributions unless the sound mode has a_pre- 
ferred wave-front orientation.’ That happens to be 
the case for progressive waves from isolated sources and 
for the diffracted pressure gradients but not for the 
other categories of sound where more random orienta- 
tion prevails. Consequently, there will remain some un- 
certainty as to the composition and intensity of free- 
stream disturbances in the tunnel when measured by a 
single hot wire no matter how accurate the electrical 
measurements may be. 

Actually, measurements of free-stream fluctuations 
in good supersonic tunnels are of the threshold type, 
plagued by problems of calibration, noise, frequency 
resolution, and unsteady strain-gage effect,’ so that the 
decomposition of the signal into the constituent flow 
disturbances is further complicated and uncertain. De- 
pending then on the assumptions of specific orientation 
of sound wave and the presence or absence of correla- 
tion between the entropy and the vorticity modes, the 
free-stream data indicated the vorticity-mode measure 
u’/ Ul, from 0.0003 to 0.0007 and the entropy-mode 
measure A7’/7\ from 0.0003 to 0.000S—i.e. 
(solenoidal) velocity fluctuations wu’ from 0.5 to 1.05 


r.11.S. 


ft./sec. and r.m.s. (nonisentropic) temperature fluctua- 
tions from 0.1 to 0.26°F. In view of the presence of 
extraneous electrical signals these results are believed to 
be high by a factor up to 3. 

Before turning to the corresponding indications of 
sound intensity, it is instructive to compare these re- 
sults with the prediction of Ribner-Tucker linear 
theory.” For the quoted wu’ Ll) of 0.0074 in the set- 
tling chamber (assumed isotropic) the theory would 
lead to the longitudinal component uw’, U’ of 0.00023 and 
the lateral component v’, LU; of 0.00028. 
measurements at low speeds (see Figs. S and 9 of refer- 
ence 27) one can infer that the nonlinear effects and the 


From Uberoi's 


viscous decay, which were neglected in the theory 
would tend to equalize the lateral and the longitudinal 
components and reduce the total energy. Thus meas- 
urements and theory, with proper interpretation, indi- 
cate that the vorticity level in the supersonic test sec- 
tion was indeed very low. As to the entropy fluctua- 
tions, one would expect these to be preserved in the pas- 
sage through the nozzle except for some decay due to 
heat conduction. Indeed, temperature fluctuations on 
the order of 0.16 to 0.42°F. 
measurements in the settling chamber (the low-speed 


were consistent with the 


temperature sensitivity being relatively poorer). The 
level of temperature (entropy) fluctuations in the 
supersonic test section was then also very low since 
fluctuations of 0.2°F., for instance, need to be compared 
with the total change of static temperature of 172°F. 
across the laminar layer. 

In view of the low absolute level of vorticity and en- 
tropy fluctuations in the supersonic section, the lack of 
shift in transition when the settling chamber turbulence 


was raised sixfold appears now less surprising. Origi- 
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nally it was inferred from the supersonic hot-wire data 


that the corresponding rise in vorticity could not have 
exceeded 15 per cent. It now appears likely that the 


assumptions of the absence of sound of categories 


and (b) and of the absence of an extraneous electrica] 


signal, which led to the 15 per cent estimate, were un 
justified. Since uncorrelated signals are added in the 
square, the low vorticity in the supersonic section may 
have risen four to sixfold and caused only a small in 
crease in the r.m.s. electrical signal in presence of th 
other larger disturbances. 

In absence of strain-gage signals the inferred sound 
intensity in the free stream ranged from 125 to 131 db 
i.e., r.m.s. pressure fluctuations were from 0.2 to 04 
per cent of free-stream pressure. The corresponding 
values of the r.m.s. longitudinal component of the ir- 
rotational velocity fluctuations “’ would be from 0.7 
to 1.4 ft. sec. if the dominant sound were of category 
When estimates of pos 


(c)—shivering Mach waves. 
sible extraneous signals are taken into account, the 
sound level may be lowered by as much as 10 db. and the 
velocity and pressure fluctuation decreased by a factor 
up to 5. 

The most accurate sound measurements were those 
of “incremental” sound—t.e., from the rise in signal 
between two points on the same streamline which is 
due to a localized increase in sound generation at the 
wall (such as that due to an element of roughness, dif- 
fracted by the turbulent boundary layer or due to 
transition on the plate). Because the background con- 
tributions are subtracted it is not necessary to assume 
specific orientations for the analysis—most often the 
near Mach line orientation is deducible from the data 
In the same way one can analyze their spectra. A 
1955 experiment by the author’ traced a very weak 
shock (1/; = 1.76, steady Ap p, = 0.049; unsteady 
pressure signal = 145 db., Ap’/p, = 0.02; total un- 
steady horizontal spread sensed by the imperfectly 
resolving hot wire = 0.004 in.) as it impinged upon the 
turbulent wall boundary layer (0.6 in. thick) and re- 
flected outward. The reflected unsteady wave ex- 
ceeded 125 db. level over a horizontal distance of 0.4 
in. and reached a maximum? of 134 db. (r.m.s. Ap’ = 
0.7 per cent of local static pressure). Its normalized 
energy spectrum as well as the vorticity and entropy 
spectra at a representative point in the reflecting bound- 
ary layer (y/6 = 0.2) are shown in Fig. 4. 

It would therefore appear that imperfections 
nozzle contours and wall giving rise to narrow 4. 
pressure changes are accompanied by a.c. pressure waves 
which spread appreciably upon each reflection from the 
turbulent boundary layer. Experience, as well as calcu 
lations (based on data of reference 38), indicates that 
most Mach waves discernible by sensitive schlieren ap- 
paratus have an a.c. content on the order of 120 to lol 
db. even if they appear steady to the eye. Many tunnels 
are then dense with shivering sound waves. One call 

+The sharp shocklet was followed by a gradual Prandtl 
Meyer expansion wave to Jf = 1.96 which may have reduced the 


intensity of the reflected wave somewhat 
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tions in a turbulent boundary layer at 14, = 1.76 and of the 


sound reflected from the laver 


speculate that their strength varies with 60(log p) /Ox, 
where p(x) would be the pressure field generated by the 
obstacle in an inviscid fluid, and with the intensity of 
the local Mach Number fluctuations. 
fluctuations would cause larger spatial excursions of 


Lower frequency 


wave fronts and would increase the sound strength. Fur- 
thermore, the various parameters governing upstream 
influence in supersonic boundary layers as analyzed by 
Lighthill*' should also affect the present phenomenon. 
Thus, for a given wall geometry this sound of category 
c) is likely to decrease with Mach Number [while that 
of category (b) may possibly increase }. 

In some tunnels stronger waves (even when weak 
enough to be missed during tunnel calibration*’) were 
shown to be the direct cause of early transition (see 
reference 14, p. 8, and reference 30, p. 5). One must 
keep in mind that lateral Mach waves from sidewalls 
of tunnels are more likely to be missed because they are 
not visible with schlieren apparatust and are sensed by 
hot wires only if mounted vertically. Transition may 
be hastened by such pressure waves through (1) local 
steady or unsteady separation of the boundary layer, 
2) inducement of steady three-dimensional boundary- 
layer flow which is less stable?* (on bodies of revolution 
generally and on flat plates in presence of lateral pres- 
sure waves), and (3) feeding of disturbance energy into 
two- or three-dimensional Tollmien-Schlichting waves. 
For the JHU boundary layer the dangerous frequency 
band was estimated to fall between 12 and 30 ke. per 
sec., and a glance at Fig. 4 shows that the relative con- 
tribution to that band from the sound of category (c) 
is on the same order as the contributions from free- 
Stream turbulence at low speeds.” '* By contrast, 


iree-stream vorticity and temperature fluctuations 
which were generated in the low-speed parts of the 
supersonic tunnel system by grids and rods*' were 
estimated to have no appreciable energy above 4-8 ke. 
and were therefore less likely to feed the 7-S waves 
in the supersonic boundary layer. 

Actually one of the many reasons for discontinuing 
the experiment was that no regular frequencies of the 


signals were discernible in 1952 and 1953 (when no 


t Their reflection from cylindrical models or probes is visible 


rion 


E 


X PERIMENTS 


suitable harmonic analyzer was available) and it was sur- 
that a different dominated the 
transition pattern. Recently, Laufer and Vrebalovich® 
discovered that Tollmien-Schlichting waves do build 


mused mechanism 


up at supersonic speeds as well, but that their ampli- 
fication rate is very much smaller than at low speeds 
(as predicted by recent advances in the stability theory), 
a fact that explains the need of a harmonic analyzer 
for their detection. The mechanism and the final cri 
for the breakdown of 


the actual 


terion 7-S waves into a tur 


bulent pattern—1.e., occurrence of transi 


remains unclear at any speeds 
a potentially strong 


tion 

One clue as to the presence of 
transition promoter came from the early flattening of 
the traverses off the flat 
plate (see Fig. 1 The measured static pressure or 
Mach number gradients, quoted earlier, did not appear 
for other for in 


surface-tube centerline on 


than those tunnels 
stance, for reference 13 where higher transition Reyn 


In particular, they were 


more serious 
olds Numbers were reached. 
weak enough not to cause steady local two-dimensional 
an axial trav- 


boundary-layer separation. However 


erse with the hot-wire vertical just outside of the 


boundary layer along 2 0.82 in. disclosed a strong 
peak in the sound signal at . = 4.4 in. Superposed on 
the gradual mean pressure field from the corner of the 
flat plate was a relatively narrow sound wave of maxi 
mum intensity of 136 db. (Ap’ pi = 0.0073, uo = 2.6 
ft. sec.). It is conceivable that this oblique wave oc 
casionally reached instantaneous values sufficient to 
separate the layer locally and cause the time-averaged 
The 


unsteady character was impressed upon the pressure 


flat variation of pitot wall pressure seen in Fig. 1. 


wave by the sidewall boundary layer so that its spec 
trum probably resembled that of the reflected sound in 
Fig. 4. The corner waves then could also feed energy 
into the unstable Tollmien-Schlichting waves. 

An attempt was made to compare the hot-wire traces 
upstream and downstream of this oblique corner wave 
but no differences were discernible with the equipment 
then available. This was indeed disappointing in view 
of the overall indications that there was some cause 
dominating the transition pattern in the JHU tunnel 
and that the corner and side-wall interference formed 
the strongest disturbance measured. This inability to 
trace the direct effect of such disturbances on the bound- 
ary layer was the final reason leading to the conclusion 
that a more versatile tunnel (allowing higher Mach 
Numbers and a wide variation in stagnation pressure) 
as well as improved equipment and techniques were 
needed for a basic experimental program on supersonic 


transition. 


CONCLUSION 


Partial description of the transition pattern on a 
plate and of free-stream disturbances in supersonic 
wind tunnels as sensed by a hot-wire anemometer has 
been presented. The specific case was used as an illus- 
tration of existing limitations on our knowledge of the 
set of conditions under which a given pattern of transi 
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tion occurs on a given model in a given testing facility. 
A particular supersonic transition experiment must 
indeed be evaluated with care and all inferences drawn 
with caution. Still, the cumulative experience from 
different facilities has led (despite many contradictions) 
to valuable empirical conjectures and to substantiation 
of some trends of transition predicted by linearized sta- 
bility theory. However, for more basic understanding 
of the transition phenomenon, broader studies (includ- 
ing features of efforts such as in references 2, 3, 17, 25, 
and 33) in the best of available facilities with sub 
stantially improved instrumentation remain essential. 
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On the Frequency Spectrum Generated by a 
Maneuvering Aircraft’ 
— 


A. P. BROGLE* 


U.S. dimy Signal Engineering Laboratories 


SUMMARY 


In this paper, an upper limit in amplitude for tne frequency 
spectrum of aircraft motion is found. The existence of a bound 


n the absolute value of acceleration——i.e., a finite centripetal 


ieceleration capability of aircraft—is the underlying phenomenon 
vhich determines the limit on the transformed position motion 
In determining this bound, only periodic motions of the aircraft 
re considered, and the only aircraft characteristic taken into 


iccount is Maximum acceleration It is shown that the derived 


maximum may be applied to the determination of the bandwidth 
direction or 


haracteristics for a communication-control fire 


interception system \ simple relationship is found to exist 
vetween the value of the amplitude limit at any frequency and the 


maximum error in position, without regard to time delay and 


noise, Which occurs in a communication-control system of speci 
fied bandwidth. Since a dominant criterion of these systems 
uppears to be to obtain adequate performance during maneuvers 
f the aircraft, the general result applies to an important opera 


tional problem. Several examples of the usefulness of the result 


in the determination of such system characteristics as band 


width, sampling interval, or transfer function for a specified 


ire illustrated 


iccuracy 


INTRODUCTION 


| Dpaennegeng and control systems in which the 
position of a moving aircraft is the applied infor- 
mation variable are receiving considerable attention 
Like 
communication 
ibility of a limited band of frequencies to convey to 


today. other communications systems, these 


fire-control systems depend upon the 


the ultimate receiver an up-to-date reproduction of the 
flight information, even though the information itself 
is not spectrum limited. This condition is, in general, 
basic to all of these systems and specifies that fine 
structure of the information is inevitably lost during 
measurement, transmission, or computation. Since 
a limited time is prescribed in which to recover flight 
information in the communication - fire-control system, 
this loss of detail always results in distortion at the 
system output. An important problem of these sys- 
tems 1s to assess the magnitude of this error in terms of 
the characteristics of the applied information and for 
given bandwidth conditions. 

With the ever-increasing problem of reduced channel 


capacity for communication-control systems, a further 
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question arises concerning the optimum balance be 
tween signal and noise to achieve best performance 
within the confined bandwidth. For example, if the 
available channel bandwidth is weighed too heavily 
in favor of noise, signal bandwidth is reduced, and the 
dynamic errors may become excessive. On the other 
hand, noise cannot be ignored to obtain rapid following 
since the fine detail obtained may be mostly noise. 

In the presence of distortion and noise then, optimi 
zation of a system results only from the use of as much 
a priori knowledge of the nature of the applied informa 
tion as can be obtained. For systems of the type ct 
interest here, the information variable is any one of an 
this 


infinite number of possible flight paths.¢ Since 


statement alone is not a sufficient description, it is 


desirable to find some common characteristic of the 
input variable that is meaningful in terms of the oper 
ational problem to be solved by the system. For 
example, if it is the average behavior of the system 
that is important, a flight path can be described, by 
statistical means, that is representative of the class ct 
input functions to be handled. This method is often 
a convenient one since noise is also best described in 
Statistical terms, and techniques are then available 
which can specify the average response for a given trans 
fer function or, alternatively, a system transfer function 
that yields best average response, usually in the mean 
square sense. However, it is believed that this and 
similar averaging methods do not, in general, satisfy 
an important operational requirement which is placed 
upon communication - fire-control systems. 
In examining the usual application of these systems 

namely, to air defense and to air navigation and con 
trol 


conditions, a primary condition sought is adequate 


it becomes apparent that, irrespective of noise 


performance during the times when evasive actions 
This 
is so because of the severe operational consequences 
usually attached to the loss of a target. Since it is at 
just the times of maneuver that the system has most 
probability of 


or turning motions of the aircraft are occurring. 


difficulty in following the target, the 
losing or of missing the target is also greatest during 
this same period. Thus, if aircraft maneuver 1s an 
event which may well defeat the purpose for which a 
system is designed and if that event is a likely occur 
rence, an essential requirement of the system 1s estab 
that a domi 


lished. It immediately follows, therefore, 


nant criterion of these systems must be to keep the posi 


t The term ‘“‘flight path” is ised to denote the continuous line 


in space that is generated by a moving aircraft 
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tion error during anticipated maneuvers of the air- 
craft within the resolution capability of the system. 
Resolution capability as used here is a rather general 
term and applies to the ability of the system to per- 
form acceptably any or all of the functions of initial 
tracking of the aircraft at the system input, of trans 
mitting this the communication 


channel, and of aiming the guns or directing the inter 


information over 


ceptor or missile to the target at the output. 


Although the foregoing is a principal argument sup- 
porting the maneuver criterion, it is well to point out 
that many factors other than those indicated are impor 
tant in specifying the detailed performance required for 
any real application. Nevertheless, the significant 
point here is that aircraft maneuver appears to be a 
most important input form to be considered for these 
systems and that adequate protection against this event 
probably has to be obtained at the expense of a higher 
level of average performance than would be required 
by a less stringent, but more unrealistic, criterion. 

If position error during maneuver is to serve as a 
useful measure of system performance, it is evidently 
desirable to find the relationship which exists between 
the position error at the system output and both the 
characteristic of the system input variable which speci- 
fies maneuver capability and the characteristics of the 
system itself. In this paper, such a relationship is 
found which appears to be useful in a number of impor- 
tant cases. 

Although several approaches to the problem of 
defining the beliavior of the applied information vari- 
able are undoubtedly available, a method which has 
possible advantages is a study of the transformed 
i.e., the frequency spectrum of a maneuvering 
such a 


motion 
aircraft. An immediate advantage is that 
description permits consistent usage of the well-known 
frequency domain language in which the communica 
tion or control system is usually described. Recalling 
that an infinite number of flight paths is possible, it is 
seen that this spectrum may take on an infinite variety 
of shapes. <A_ possible limitation of the frequency 
domain technique in expressing a priori information in 
concise form therefore becomes apparent, especially 
since statistical techniques which lead to average results 
do not appear to best match the physical situation 
here. However, it will be shown that a particular 
spectrum can be found which defines the upper limit 
in amplitude under which all spectra of moving air- 
craft of specified acceleration capability are constrained 
to fall. Furthermore, this spectrum is obtained for a 
maneuver that appears to be a highly probable form of 
evasive action. Thus, a common characteristic of the 
input variable which specifies maneuver capability is 
determined. It will also be shown that the amplitude 
limit at any frequency is related simply to the maximum 
“bandwidth” error of the communication or control 
system which attempts to follow the continuously 
changing position of a moving aircraft. A powerful 
link is therefore established between a fundamental 
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characteristic of aircraft motion and communication 
control system performance. 


AMPLITUDE LIMIT 


In general, the Fourier series furnishes a method 6; 


obtaining an analytical expression for a function be 


tween any two arbitrary values of the independent 


variable. 
the path of a moving aircraft over any arbitrary dura 
tion of time.* In this manner, the position in spac 
of a moving aircraft at any time / is represented by 
k2(t) = 


r(t) x(t) + jy(t) + 


/ 


i >> X(nje"™ +7 D> YV(man)e’™ 


k > 2Z(new)e’” 


where 7, 7, k are the unit direction vectors in Cartesia 
coordinates. X (nw), Y(nwo), Z(mwo) in this expressior 
are the coordinate frequency spectra of the moving air 


craft and are related to the original flight path r(t) by 


>] 
X(nwo) = (1/7) x(the 7"*" dt 


27 
Y(mwo) = (1/7) / vibe gd? Z 


« 


+7 
Z(nay) = (1/7) sie" di 


0 


* Although the Fourier integral permits a more conveniet 
representation in some respects, the infinitesimal amplitudes 
the transformed motion in this method present certain difficulti¢ 


in interpretation not encountered in the series representatio 


where amplitude is preserved 
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where / 2n/wy is the arbitrary expansion interval 
wer Which the aircraft motion is of interest and out- 
side of which the motion is assumed to be repeated. 
Differentiation of Eq. (1) gives the velocity and acceler- 


ition of the aircraft as 


) X(t Jv(t) + R(t) 
S ( jnwy) X (nwo)e?"™" + 7 D> (jnwo) ¥(nax)e?"™ 4 
k > (jne)Z(nex)e?"" (3) 
r(t) v(t) + H(t) + R3(t) 
> — (nay)? X (nape?"" + 7 
nw)? ¥(nu)e2"™ + k 2 — (na)?Z (nowy) " (4) 


The coordinate transforms of Eq. (4) are 


X, (nwo) — (nay)? X (now) = 
>! 
ci f3 | (pe “"" dt 
Y, (nan) = — (nao)? V¥(nan) = 
“7 fons (5) 
(1/7) | p(the 7" dt 
Z, (nan) = — (nwo)?Z(nwo) = 


+1 
(1/7) B(t)e 7" dt 
By revising Eq. (5) slightly, the spectrum of r(t) may 


now be written as 


+] 

X (na) = [1 /(nw)?T #()e~7"™ dt 
a7 

V'(nw) [1 (max)? 7°] | yHe 7" dt (6) 
+1 

Z(nw) = — [1 (nus)? | a(tje uate 


These equations express the frequency spectrum 
generated by a moving aircraft as a function of the 
coordinate accelerations X(t), (t), and 3(¢t). These 
components are referenced to the coordinate system 
from which the motion of the aircraft is observed, 
rather than the coordinates of the aircraft where the 
accelerations causing maneuver are actually applied, 
and it is desirable to relate the spectra of Eq. (6) to 
the vector acceleration 7(t). This is done quite simply 
by noting that the vector acceleration at any instant 
of time is defined by its projections on the fixed coor- 


dinate axes: 


x(t) = rit) cos @ 
y(t) = r(t) cos B (7) 
g(t) = r(t) cos ,| 


where a, 8, and y are the angles formed between the 
vector acceleration and the x, y, and z axes, respectively. 
hus, if Eq. (6) is considered as defining a vector 


spectrum 


EK 


i es 
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R(ne) = 1X(noo) + 7) kZ (na) 8) 


Mw) T 


the transformed motion, as a result of Eq. (7), becomes 


R(nw) = —[1/(new)?7] X 


2] e/ 
E r(t) COS ae "at+ J r(t) COS pe dt 
* . 


The frequency spectrum R(jia) is now stated in terms 
of the acceleration resulting directly from the forces 
produced on the aircraft For convenience of termi 
nology, this acceleration will hereafter be called the 
“applied” acceleration. 

fundamental statement 


For aircraft, however, a 


concerning /(f) is that its values are bounded—-namely, 


r(t)| <A (10) 


where A is the maximum acceleration capability of the 
aircraft. It is important to note that this condition 
specifies only a maximum value for 7(¢) and does not 
restrict the variation of vector acceleration in any way 
within this range.* The 
determine the manner in which the restriction on the 
range of r(¢) affects the magnitude of the position spec- 


main problem here is to 


trum of Eq. (9). 
The result of R(na ) because of the existence of a 
finite bound on the magnitude of the vector accelera 


tion can be shown to be 


R(nw)mar = } ~4A/m(nwn)?; X 
[1 cos a + j cos 8 + kReos y] (11) 
where » now takes on only positive values. It can 


also be shown that maximum R(nwo) is achieved only 
when /(/) is an applied square wave of acceleration of 
amplitude A and period 27/(mw»). In practice, this 
maneuver is obtained in an aircraft by the repeated 
application of maximum acceleration first in one direc- 
tion for a time 7/nw) and then in the opposite direction 
for the same length of time. The resulting motion is 
a more or less weaving or serpentine path consisting of 
a series of equal left and right turns of the aircraft or 
a prolonged climbing and diving path (as illustrated 
by the curves of Fig. 3). This square wave of acceler 
ation both satisfies the condition for complete variation 
within +A set down previously and appears to approxi- 
mate a real case of evasive maneuver. The range of 
Eq. (11), therefore, not only is a theoretical maximum 
but can serve as a practical maximum as well. 

In order to visualize the conditions of Eq. (11), it is 
helpful to observe that the terms within the square 
brackets define a direction vector of unit magnitude 
corresponding to the instantaneous direction of the 
These terms relate the orienta- 
the fixed 
j, and k. 


applied acceleration. 
the coordinates of aircraft to the 


The term out- 


tion of 
cocrdinates described by 2, 
* The rate of change of acceleration is also limited because of 
aerodynamic considerations. The corrections to the basic analy 
sis due to the finite rate of an aircraft's response to the controls 


are indicated in the Appendix 
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Fic. 2. Applied maneuver to attain bound of Fig. 1 


side the square brackets represents the maximum ampl1- 
tude that can be achieved by any sinusoidal frequency 
Taken 


together, these two terms define a vector frequency 


during a maximum maneuver by the aircraft. 


spectrum whose maximum value in any of the fixed 
coordinate directions is continuously modified by the 
projection of the acceleration-direction vector on that 
coordinate. However, a simpler interpretation is 
possible when only the magnitude of the position spec- 
trum is considered. 

Since the absolute value of the bracketed terms of 
Eq. (11) is by definition unity, the magnitude of 


R(nw,) becomes 
R(na) | mar = 4A/a(nwo)? (12) 


Therefore, with reference to Eq. (9), the maximum 
amplitude of any frequency component of the posi- 
tion variable r(/) can never exceed the value expressed 
by Eq. (12). Although this amplitude limit can be 
treated as a scalar quantity, it is important to remem- 
ber that the limit at any frequency is actually free to 
occur in any direction. The spectrum of Eq. (12) is 
thus a spatial representation of the upper bound of the 
frequency spectrum generated by a maneuvering air- 
craft. 

The importance of Eq. (12) in limiting the frequency 
spectrum of aircraft motion is demonstrated by the 
curves of Fig. 1. In the Figure, it can be observed 
that frequency components that are in the vicinity of 
0.1 eps, for example, are permitted a maximum ampl- 
tude of approximately 65 yards for a 2g acceleration, 
whereas frequency components near 1.0 cps are allowed 
an amplitude of but 0.65 yard for the same limitation 
on the acceleration. Therefore, the component ampli- 
tudes of motion in space are significant only over a 
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rather small range of the frequency variable (ne, 
The extent of this band which must be preserved jp 
order to meet given system error requirements will be 
the subject matter of the remainder of this paper. 


RELATION OF FREQUENCY SPECTRUM TO POSITION ERROR 


In the previous section it was stated that an applied 
square wave of acceleration is the only form of motion 
which produces amplitudes equal to the bound of Fig 
1. Fig. 2 shows the essential characteristics of the 
applied acceleration motion and the corresponding 
velocity and position variations resulting from this 
disturbance. As before, the square wave of accelera 
tion is taken only with respect to the motion of the 
aircraft. The solid curves represent the total motion 
while the dashed curves indicate, for illustration pur- 
poses, the fundamental and third harmonic compo 
nents. Also shown are the corresponding analytical 
expressions for the motion. Integration of the applied 
acceleration produces a linear variation of velocity 
and integration of the velocity yields the resultant 
parabolic position variation. Since the main interest 
here is in the motion due to maneuver, the constants 
of the integrations are not included. The frequency 
components produced by the successive integrations 
occur with the amplitudes and phases shown in order 
that their sum always produces the total motion 
With reference to the resultant position variation, it is 
important to note that, because of the rapid conver 
gence of the position series, the fundamental component 
of amplitude A / 7*f*,* is not significantly different per 

* For the remainder of this paper the continuous variable 


can replace the discrete variable (mw) without loss of generality 
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centagewise from the total motion of amplitude A 52/ 
In fact, it is only at rather low frequencies that the 


difference in terms of the absolute distance becomes 


ippreciable. 


In order to illustrate the foregoing statements and to 
lemonstrate by example the validity of the spatial 
representation of the previous section, a two-dimen 


aircraft motion is considered. Fig. 5 


4% possible flight paths occurring 1n a 


sional case of 


shi WS 


i number 


horizontal (y-y) plane for various values of the fre 


of the applied square wave of acceleration of 


quency 


my ) 
1g 


4 _- 


in the general direction of the x coordinate at different 


The uppermost four curves represent flight 


lrequencies of maneuver, while the lowest curve is an 
example of a flight path with an average direction 45 


to the v and y axes and for the same frequency as the 


fight path immediately above it. A maximum accel 


eration of A 2 and 


—S 


a velocity of 500 m.p.h. have 


been chosen for the illustration. Since the motion in 


the x-y plane will be circular, the corresponding values 
4 turn radius and turn rate for the paths are RK = 5,000 


ids and w = 1S rad. sec., respectively. Time 1s 


shown by the dots which are placed at 5-sec. intervals 


long the curves. Similarly, the vector acceleration 


is shown for a portion of the lowest curve by the 


irrows, all of which are of magnitude A. Figs. 4 and 


» show the resulting coordinate acceleration and posi- 


iriations as functions of time. If a careful 


xamnation is made of the curves, it is seen that, even 
the motion is relatively complicated in its 
mean 


letails, the amplitude of the motion about the 


position 1s almost exactly equal to the corresponding 


nplitude given by Fig. 1 for all frequencies above 
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approximately 0.025 cps.* For frequencies below 


0.025 eps and with reference to any particular coord: 
nate, the amplitudes of Fig. | appear to be high by a fac 


tor of approximately 43. This condition 1s as expected, 
however, since the applied coordinate acceleration at 
these frequencies is nearer sinusoidal than square wave 
in shape, and Fig. 2 shows that the resulting position 
motion for a purely sinusoidal acceleration will have 


an amplitude 7 4 as the amplitude for the 


as large 
square-wave case. Since this apparent difference 0¢ 
curs well within the pass-band in most applications 
and therefore lies outside the range of interest here, 
and because it does not affect the spatial representation, 
it is not included as a refinement to the curves of Fig. | 

Examination of the possible conditions existing at a 


particular frequency / of Fig. | also permits a valuable 


observation to be made. At the outset, it 1s evident 
that a component can occur at f due either to the 
presence of a fundamental component resulting from 


at 


the application of the square wave of acceleration 
that frequency or to the presence of a harmonic of an 


equivalent motion whose fundamental is at some lowe! 


frequency. If the fundamental 1s present at /, its 
maximum amplitude can be Az On the other hand 
the maximum amplitude that a third harmonic can 


A/2ix 3 
the maximum amplitude of a 


have at this frequency is, from Fig. 2, just 
or A/3r*} Similarly 
fifth harmonic 1s A 57°*f*, and further harmonies take 


the approximate form A nq Since the foregoing 
irgument holds at every frequency, it is seen that the 
fundamental component of the motion alone determines 
the bound of Fig. 1. 

Phis frequeney will differ slightly for different velociti 


However, it can be verified that velocity is not a significant f 
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It is now appropriate to pursue the main problem of 
relating the special properties of the frequency spec- 
trum of a maneuvering aircraft found by the previous 
analysis to the performance of a communication or 
control system. Stated simply, the basic operations 
performed by such a system are usually (1) to measure 
and track the path generated by a moving aircraft, (2) 
to recover the path at a remote location or to perform 
some other operation on the information, and (5) to 
carry out these operations accurately, in real time, and 
without undue time delay. A fundamental character 
istic of any of these systems is the availability of a 
limited band of frequencies in which to measure, com 
pute, or convey flight information. 

Consider Fig. 6 which is an elementary block diagram 
of a system incorporating features of both communica 
tion and automatic control. This Figure represents 
one of a number of ways in which aircraft ““B’’ might 
be automatically controlled from the ground to inter 
cept a target, aircraft “A.” It is assumed that the 
computing facilities are located at some distance from 
the tracking site, that high precision is required, and 
that the interconnecting link may be noisy. The 
use of coded transmission is therefore probably called 
for, and sampling makes it possible to time-share the 
position coordinates of both aircraft over the same 
narrow-band channel. The computer compares the 
positions of target and interceptor aircraft and cal 
This 


signal is then transmitted by radio to the intercepting 


culates a command signal for the interceptor. 


aircraft in such a form that the automatic pilot com 
pletely takes over the flight. 

The response of a linear system, such as Fig. 6, to 
a position input r(f) will be the product of the frequency 
spectrum generated by the moving aircraft and the 
transfer-function spectrum of the system. Since there 
will always be some frequency f. beyond which the 
system does not sensibly respond to the input, a por 
tion of the frequency spectrum R(mw)) of Eq. (8) can 
be lost. In particular, if R(#w») occurs entirely above 
f., essentially no motion is observed by the system. 
Again, it is the frequency components of the maneuver 
that are important. It is reasonable then to postulate 
that the maximum error of the system, exclusive of 


noise and time delay considerations, in describing the 
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position of the moving aircraft occurs when the applied 


maneuver of Fig. 2 takes place just outside /..* J; 
this case, the maximum error in position 1s equal to thy 
total motion of f,. J 


error for this condition is 
Bae = A/se 


Since this error is not essentially different from thy 
maximum amplitude of the fundamental component of 
the motion at fo, Fig. 1 can now serve not only as the 
upper bound of the frequency spectrum generated by 
the aircraft but also as the error vs. bandwidth measur 


of the system that accepts the position variation of this 


aircraft as input. - 

Eq. (13), although derived heuristically and fron 
essentially a steady-state analysis, nevertheless pr 
vides a fundamental relationship between the physical 
behavior of a moving aircraft and a basic characteristi 
of the associated communication or control system 
Since the motions which determine the bound of Fig 
also seem to be highly probable maneuvers of the air 
craft at the time when the system is called upon to per- 
form most accurately, a realistic figure of merit i 


if 


evaluating performance is introduced. [1 


addition, it appears that the techniques which lead t 


system 


this result can be extended to take account of time 
delay and noise and can be applied to the solution oi 


It has been verified by calculation that this error is great 
than the error which can occur in a typical low-pass system wh« 


the fundamental component of the motion of Fig. 2 occurs bel 
‘and only harmonies of the motion are eliminated by the svster 
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From Fig. 2, the resultant spatial 
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other problems involving aircraft maneuvers. 


fnal section will be devoted to two examples of the 
usefulness of the results in the determination of impor 


tant communication or control system characteristics. 


DETERMINATION OF SYSTEM CHARACTERISTICS 

The svstem transfer function of Fig. 6 is a rather com 
plicated expression consisting of the closed loop rela 
tionship of the individual transfer functions of the track- 
ing units, the sampling and coding processes, the decod 
ing the computer operation, the transmis- 


sion links, and the autopilot-interceptor combination. 


processes 


On a system basis, the design of each of these parts must 
be related intimately to the performance of the others 
ind to the characteristics of the applied information. 
In general, only one of these components will be instru 
mental to specifying the cutoif frequency /. for the sys- 
tem. Application of Eq. (15) will then determine the 
maxilunum position error in space expected for the cutoff 
and the maximum acceleration capability 
On the other hand, when more 


frequency 
issigned for the aircraft. 
freedom of design is permitted, /. can be adjusted in 
order to meet the error limit set for the system. Thus, 
Eq. (13 


tem performance. 


can be used directly to predict overall sys 
A specific example is the calcula 
tion of the performance of the samplers of Fig. 6. 

The operation of sampling is usually employed in 
a communication or control system when the data are 
inherently obtained in this form or when it 1s desired 
to time-share signals from two or more separate sources 
over the same channel or to use a digital computer in 
the system. In order to preserve the information, it 
is well known that samples of the time-varying input 
must be taken at a rate that is at least twice the highest 


significant frequency of the input or 


> 2f (14) 

In terms of a uniform sampling interval 7,, Eq. (14) 
becomes 

/ * 1/2; (15) 


It is evident immediately that the highest significant 
frequency f, in the case of a moving aircraft is just the 
cutoff frequency defined previously. Therefore, from 
Eqs. (13) and (15), the maximum error associated with 
a particular sampling interval 7’, is 


Ema Al//% (16 


his equation is plotted in Fig. 7 for several values of 
the acceleration capability A. When the sampling 
operation is the bandwidth determining factor, the 
curves of this Figure express directly the inherent pen- 
alty attached to the choice of a particuiar sampling 
interval in the system. Note that an added penalty ex- 
ists for the use of longer intervals because of the squared 
nature of the variation. In cases, such as Fig. 6, where 
the sampling operation is probably not the bandwidth 
determining factor but is included as a component 
transfer function of the system, the same general tech- 


The 


E 


Ni 


SPECTRUM $4.3 


nique applies. Here, however, the sampling interval 


is chosen by Eq. (16) or Fig. 7 so that undue error is 


not contributed by the operation. In systems where 
the information is encoded into digital form either for 
transmission or for computation, these relations are 
also helpful in balancing the quantizing error (noise 
against the error or position uncertainty due to ma 
neuver. Since quantizing error decreases with longer 
sampling intervals because of the inherent smoothing 
maneuver error increases rather 


obtained whereas 


rapidly, the choice of a compromise interval is often 
quite important. 

Another example of the usefulness of the concepts 
presented in this paper is the calculation of the shape of 
the system transfer function which provides best pro 
tection against noise without degrading the accuracy 
for maneuvers beyond the limit set by E,,,,. It was 
pointed out earlier that maximum position error occurs 
whose fundamental variation 


the 


only for a motion of Fig. 2 
the cutoff 
for other conditions of motion below 


is located at Furthermore, errors 


are quite small 


compared to the maximum motion. This is largely 


due to the fact that the fundamental, and often 
several harmonics, are contained within the transmis 
sion band of the system for motion below If now 
the maximum error can be tolerated for any fre 


queney of maneuver of the aircraft, it is possible to 
minimize the magnitude of the transfer function of the 
system everywhere within the transmission band to 
a value compatible with this error. By providing min 
imum gain for the system, the transfer function is, in 
an intuitive sense, automatically adjusted to furnish 
best protection against noise. 

Let E,,,; be the maximum error in distance that can 
t) the 
In order to appor 


be tolerated for a particular system and E( per 
centage error at any frequency 
tion the error over the transmission range, the percent 


age error must be 
E(f) RAG 3 ¢ 6 (A /32} A/xf?) = 


iz iS, aS before, the maximum amplitude 


R(f) 


of motion at f—i.e. 


where 
, the upper bound of Fig. 1 at that 
frequency. If 


IS) 


is the voltage transfer function of the system, the per 
centage error at any frequency is related to this quan 
tity by 


H(f) 


(19) 


Substitution of Eq. (17) in Eq. (19) immediately gives 


H(f) = 1 — Enmaz/|R(f) | mar = 1 ff 20) 
This equation then is the desired system transfer func 
tion to allow Ener anywhere in the transmission range 
and hence to achieve noise protection. Fig. 8 shows the 
ilready existing at 


* This equation neglects the small error 
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Fic. 8. Optimum transfer function and resulting error curve to 
achieve maximum noise protection in system 


resulting system characteristics in graphical form. 
Although the transfer function of the Figure is not 
realizable exactly, it can be approximated satisfactorily 
by a few network components. It can also be observed 
that the effective attenuation of noise is greatest for 
higher frequencies in the band. In fact, at 0.9/, it is 
14 db. and increases theoretically to infinity at /,. 


APPENDIX 


Another factor limiting the upper bound of the 
frequency spectrum of aircraft motion is the restric- 
tion of a finite rate of change of acceleration correspond 
ing to the existence of control surfaces and aerodynamics. 
Thus, it is seen that the square wave of acceleration of 
Fig. 2 is not a completely accurate representation of the 
physical situation. Although the transfer function 
expressing the rate of the aircraft's response to the con- 
trols is complicated, involving many time delay fac 
tors, this expression can be approximated for illustra 
tion purposes* by the single time constant term 


G(S) = (0p/6)(S) = 1/01 + TS) (21) 
where 
6p = direction of motion of the aircraft 
6 = control-surface deflectionTt 
7 = lumped time constant of the aircraft response 
S = complex frequency variable 


Since #p, the change in direction of motion or amount 
of turn of the aircraft, is proportional to the magnitude 
of the centripetal acceleration causing turn, Eq. (21) 


also expresses the rate of change of acceleration in 


* Over restricted ranges of frequency, representation of the 
acceleration spectrum as related to the characteristics of the 
urcraft by a second-order transfer function is possible but 
literally never by the first-order function chosen here 

+ The time lag introduced by the controls themselves are 
issumed to be negligible compared to the lags of the aircraft 


response 
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Fic.9. Correction factor to Fig. 1 to take account of aero 
considerations 
response to control-surface deflection. This th 
the form of the factor which must be applied to the 
amplitudes of Fig. | of the text. 
Combining this factor with Eq. (12), the resultant 


expression for the maximum amplitude is 
>] 
R(w) = | 1A W\W)-j* | | \ | T 7 Ww 2) 


For aircraft, 7’ appears to take on a fairly wide rang 
of values depending on the characteristics of the par 
ticular aircraft, including the peak acceleration cap 
bility. Fig. 9 is a plot of the correction factor term 
of Eq. (22) for values of the response time constant of 
From the Figure, it can be noted 


1, 2, 5, and 10 sec. 
that the effect of aerodynamic lag is, as expected, mos 
pronounced at the higher frequencies and produces 
significant reduction in amplitude for frequencies in the 
range of 0.1 to 1 eps and for the higher values of 
Thus, an adjustment may have to be made in the curves 
of Figs. 1 and 7, but not in the basic method of the 
report, in applications where this difference will causé 
a marked change in the results and when the aerod\ 


namic response characteristics of the aircraft are 


known. 
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\ LJ | : P : 0 = conditions at reference station x = a 
\ | | [wo similarity laws are known for the mean-velocity profile in | dacnll aay 7 
|| oar ; : } Re = conditions at edge of boundary layer 
j turbulent boundary layer with constant pressure These are titles ok ei 
Prandtl’s law of the wall and Karman’s momentum-defect law 
\ |] he first law has recently been generalized empirically to flows 
0A i0| with arbitrary pressure gradient by Ludwieg and Tillmann, and (1) STATEMENT OF THE PROBLEM 
rod the second law to a certain class of equilibrium flows by F. Clau =) Retvediension 
In the present paper it is shown that the pressure distribution OR FIFTY YEARS the problem of turbulent shear flow 
= corresponding to a given equilibrium flow can be computed by has occupied a prominent place in the literature of 
: =: ssuming a certal ari ‘ter D = ))dq/dr, is consti . ‘ : ‘ : 2 : : 
d to tl wren he hqnntans penenmener © 4: POET Se Seen, fluid mechanics. During this time, the principle of 
i where g and r, are the dynamic pressure in the free stream and 


») 


ide range 


the par 


the shearing stress at the wall, respectively. The hypothesis 
D = constant is suggested by a study of the integrated continuity 
the 
class of equilibrium flows defined by Clauser and the class of 
Falkner and Skan 


constant is also verified using experimental data for several 


equation and is supported by a rigorous analogy between 


laminar flows studied by The hypothesis 


D = 
equilibrium turbulent flows and is interpreted physically from a 


kinematic point of view 


physical similarity has always been an important tool 
in the study of easily observable quantities, such as the 
mean-velocity distribution. For the turbulent bound 
ary layer, in particular, appreciable progress has been 
made in dealing with what I will call the direct prob- 
lem—that of describing empirically the mean-velocity 


profile in flows which have been observed experi 


; is Two hypothetical limiting cases of equilibrium flow are de mentally. However, no really satisfactory solution 
M1 Capa " ‘ ; i . , . ‘ 
tor ter aimee. At con cnt » ie Geeiey Tapet Ses ee Bw, has yet been found for what I will call the inverse 
) erm } . es io ' oc} -ofile ide » ‘ 
with a completely logarithmic mean-velocity profile outside th > — . : 
ate ew oe ee ee. re ee problem—that of predicting the behavior of turbulent 
stant ol sublayer \t the other extreme is a continuously separating , : é 
' ' , a . . ma flows which have not necessarily been observed but 
Me noted boundary layer in a dimensionless pressure gradient (x/g)dq/dx - 


ed, most 
duces a 
es in the 


] 
eS a Z 


ipproximately twice that for the corresponding laminar flow 
[vpical shearing-stress profiles are computed for several equi 


librium turbulent flows, including the two limiting cases 


SYMBOLS 


which develop in a specified environment. 

From the beginning, most attempts to solve the 
inverse problem for the turbulent boundary ‘aver have 
relied on approximate methods which are known to be 
successful in the corresponding laminar problem or in 


€ curves Naa “3 
Ps ee oe! ve the problem of free turbulence. These methods 
1 of the D = parameter defined by Eq. (5 , nga 
1 caus } = exponent in m,x7" = constant assume that the mean motion can be realistically 
eal = static pressure described by one or more local similarity laws involving 
ACT OC " . > ° 7 ° ° ae 
yep = parameter defined by Eq. (6 certain characteristic parameters. The nature and 
‘ are mee ee a Siete - ; 
pui’/2 = dynamic pressure number of the laws and parameters in any particular 
= velocity defined by Eq. (47 : ss 
oe : case depend on the data of the problem. For example, 
= mean velocities tangential and normal to wall } ; ; ‘ 
re . in a turbulent free shear flow the mean motion is 
=V°r p = friction velocity 7 " . : : 
= coordinates tangential and normal to wall usually assumed to be independent of the fluid viscosity, 
aaa = ku;/u;, except following Eq. (1 in which case the width of the shear flow and the 
= boundary-layer thickness velocities at the boundaries or at points of symmetry 
6* = displac . icknecc finer ; 299 t se o 
nena: enemanes Saeenne Sy iy. {a are the natural parameters of the flow. Similarity 
= «ku/u, or y/6 according to context = ; : ; 
e an : —ro laws for mean velocity and turbulent shearing stress 
6 = momentum thickness defined by Eq. (23 , - pide aa 
= Karman’s constant in Eq. (1 may first be derived by an elementary dimensional 
aati A = —u,dx/du, analysis. The laws governing the development of the 
- u = viscosity profile parameters with distance then follow from the 
v York = Kinematic viscosity equations of mean motion in differential or in integral 


= profile parameter in Eq. (1 
= density 
T = shearing stress 


and (25 


= functions of x defined by Eqs. (24 


= functions of 7 and y/6é defined in reference 4 
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form. A summary of the main conclusions of this 
method for several free turbulent flows has been given 
by Squire.! 

If these same methods are applied to a shear flow 
bounded by a wall, however, a difficulty arises immedi- 
ately. The influence of viscosity in the sublayer intro- 
duces a new characteristic length, say v ™, and the 
finite shearing stress at the wall introduces a new 
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characteristic velocity, say u,. But with three lengths 
y, 6, and vu, and three velocities u, m,, and u, to be 
accounted for, the most that can be said a priori is that 
uu; expressed as a function of y 6 may depend also on 
uu, and 6u,/v. Thus it is necessary to solve the 
direct problem—i.e., to describe the velocity profile in 
much more specific language than that of dimensional 
analysis—before integral laws for the conservation of 
mass, momentum, or energy can be applied with 
confidence. 

If the flow is laminar, the methods of dimensional 
analysis can often be replaced by a genuine similarity 
argument based on the equations of motion. In 
particular, similarity laws for velocity and shearing 
stress need not be formulated separately, as either 
determines the other through the Newtonian stress law 
7 = wOu/Oy. The real justification for the assumption 
of similarity in laminar flow, however, is not often 
emphasized. It is the empirical fact that similar flows 
are apparently the ones which tend to occur in nature 
whenever the physical conditions are appropriate. 
Obvious examples, in addition to various free shear 
flows, are the Blasius flow, the plane stagnation-point 
flow, the laminar free-convection layer on a vertical 
wall, the asymptotic case of large suction, and so on. 
Hence the common assumption that laminar boundary- 
layer profiles in general form a single-parameter family, 
in which w/a, as a function of y/6 depends only on the 
ratio of 7, to uu; 6, or in the.present notation of mu, 
to 6u,/v. Not only is this assumption a natural one 
but its limitations are intuitively clear once the physical 
processes involved in the generation, transport, and 
diffusion of vorticity are understood. 

Any discussion of turbulent flow, on the other hand, 
must proceed in the vacuum created by lack of under- 
standing of turbulent transport processes in even the 
simplest turbulent motions. In these circumstances 
there is usually little choice except to abandon any 
formal analogy with laminar flow and to rely openly on 
empirical similarity laws. The similarity laws in 
question have evolved very slowly, and the efforts of 
many investigators have been required to reach even 
the present imperfect state of the art. I think it is 
important to remember, however, that the similarity 
laws to be cited here were originally discovered by a 
study of turbulent flows which also occur, so to speak, 
in nature—namely, the fully developed pipe and 
channel flows and the turbulent boundary layer with 
constant pressure. 

Unfortunately, it is usually found in practice that a 
solution of the direct problem for the turbulent bound- 
ary layer, although necessary for the application of 
integral methods, is not sufficient to allow a solution 
of the inverse problem. In fact, the only analytical 
solution worthy of the name which has so far been 
found is Karman’s analysis of flow at constant pressure. 
The difficulty is simply that any more complicated 
flow than the one treated by Karman must involve at 
least one additional parameter in the mean-velocity 


profile. At whatever level a solution is attempted, 
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some theoretical or phenomenological argument has t 
be found to connect this parameter with other quant; 
ties describing the flow. The search for such a principk 
has been carried out in many directions, and accounts 
for a large part of the extensive literature in this field 
including the present paper. 

No attempt will be made here to solve the inverse 
problem for an arbitrary turbulent boundary layer. | 
propose instead to study certain turbulent flows which 
are analogous to the laminar flows of Falkner and 
Skan in the same sense that the turbulent boundary 
layer with constant pressure is analogous to the laminar 
flow of Blasius. The main element in the analogy is a 
new concept of kinematic similarity—1.e., of similarity 
based on the continuity equation. This concept will 
be introduced following a brief but explicit statement 


of the inverse problem. 


b) The Direct Problem 


Various theoretical and experimental investigations 
have led long ago to the recognition of two similarity 
laws for the mean-velocity profile in the turbulent 
These are Prandtl’s law of the wall 
Recently both 


boundary layer. 
and Karman’s momentum-defect law. 
concepts have been extended empirically, the law of 
the wall to flows with arbitrary pressure gradient by 
Ludwieg and Tillmann,’ and the defect law to a special 
class of equilibrium flows by F. Clauser.* 

Taking the work of Clauser as a point of departure 
I have proposed in another paper‘ a formulation of the 
mean-velocity profile which incorporates the similarity 
laws in their extended form but which is not restricted 
to equilibrium flow. To be specific, consider turbulent 
flow of an incompressible fluid past a smooth plane 
surface, at which the velocity vanishes and the friction 
is Newtonian. For flows which are steady and two- 
dimensional in the mean, it is found that the mean- 
velocity profile u(x, y) may be quite generally repre- 


sented by a formula 


u(x; ¥) | yur(xc) | . w(x) y 
=f es wi ] 
u,(x) i K 0(x) 
where u,(x) = V ry(x)/p. The quantity 7,,(%) is the 


wall shearing stress; 6(x) is a boundary-layer thickness 
which is uniquely defined; and « is Karman’s universal 
constant, taken here as 0.400. The function f(yu, ¥ 
called the law of the wall, has the properties that f(z) > 
s for s > 0 (<1, say), and f(z) — (1/x) Inz + c, where 
c is a second constant taken here as 5.10, for zs > © 
(2>50, say). The function w(y ‘6), which I have called 
the law of the wake, has by definition the properties 
that w(0) = w’(0) = 0, w(1) = 2, and f, (y d)dw = 1. 
The parameter 7(x) is then a measure of the relative 
magnitude of the wake and wall components of the flow. 

I want to emphasize that the two functions called 
the law of the wall and the law of the wake are to be 
regarded as completely empirical functions established 
by direct observation of the mean-velocity profile. No 
attempt is made here to discuss the problem of turbu- 
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lence per se, and this omission is both the greatest 
and the greatest 
development. At the same time, the analogy which is 
implied by the term “law of the wake’ is meant to be 
taken seriously, in that a knowledge of turbulent mixing 
processes in free shear flows should have direct applica- 


tion to the more complicated problem of the turbulent 


strength weakness of the present 


boundary layer. 

According to the experimental evidence, * 
tions under which the two-parameter formula (1) is 
valid are nearly equivalent to the conditions under 


the condi- 


which the boundary-layer approximation itself might 
be expected to apply. The formula will certainly fail 
in the vicinity of strong anomalies in either the surface 
geometry or the boundary conditions—for example, 
close toa spoiler. On the other hand, there is no reason 
as yet to suppose that the formula must necessarily fail 
in flows which are actually separated from the adjacent 
wall. The important point is that Eq. (1) provides a 
quite accurate analytical representation of the mean- 
velocity profile for a large class of flows as a linear 
combination of two universal functions f(yu, vy) and 
wy 6). Eq. (1) therefore constitutes a useful if 
tentative solution of the direct problem for a general 


turbulent boundary layer. 


c) Equilibrium Flow 


A local friction law for the turbulent boundary layer 


is obtained from Eq. (1) by putting y = 4, 


uu. = (1/x) In (6u,/v) + ¢ + (27 /k) (2) 


and the defect law is then obtained by forming the 
difference 


(wu) — u) u, = —(1 x) In(y 6) + (@ «)[2 — w(y/d)] (3) 


or, more generally, 
(u, — u)/u, = F(a, y/b (4) 


If r is constant in Eq. (4)—1.e., independent of x—the 
momentum defect (uw, — wu) u, depends only on y/6 
outside the sublayer. This property defines a turbu- 
lent equilibrium boundary layer for the purposes of the 
present paper. It may sometimes be convenient to 
exhibit the fact that the defect law (4) varies from one 
equilibrium flow to another by retaining the notation of 
Eqs. (1) to (3), in which the parameter 7 is defined 
explicitly, without, however, assuming that the wake 
function w(y 6) has any universal validity. In any 
event, the defect law will be found to be an essential 
element in the later discussion, while the law of the 


wake will not. 


d) The Inverse Problem 


If the similarity laws of the previous sections are 
valid, the dimensionless mean-velocity profiles form a 
two-parameter family whenever the flow is bounded 
by a wall. This property supposedly expresses the 
simultaneous importance of turbulent and laminar 
transport processes in different regions of the flow. As 
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pointed out in the Introduction, problems in free 


turbulence can be viewed as a special case for which the 
parameter wu, and also the kinematic viscosity vy may 
be omitted. The number of dimensionless parameters 
is then reduced, the integral conservation laws are 
easily applied, and various well-known properties of 
free turbulent shear flows—e.g., the linear rate of 
growth of a turbulent jet 


The development of a general turbulent boundary 


are easily inferred. 


layer, however, can only be completely described in 
terms of two constants uw and p characterizing the fluid 
and four parameters , “,, 6, and m characterizing the 
state of the flow, the latter quantities being considered 
as functions of a single variable x. Knowledge of these 
four parameters in any region determines not only the 
surface friction and the rate of boundary-layer growth, 
but the complete mean-velocity profile and, therefore, 
at least within the limitations of the similarity laws 
and of the usual boundary-layer approximation, the 
shearing-stress field and the rate of energy transfer 
from the mean motion to the turbulent secondary 
motion. 

Four independent 
parameters ™, u,, 6, and mw are obviously required in 


relationships among the four 
any general formulation of the inverse problem founded 
on Eq. (1). Two of these relationships are provided 
by the local friction law (2) and by the Karman mo 
mentum-integral equation which will be introduced in 
its proper place. A third relationship is ordinarily 
included in specified conditions for a particular flow 
e.g., 7 = constant for an equilibrium flow, or #, = (x 
for a prescribed external flow. 

The assumption that a fourth equation can be found 
is not essentially different from the traditional single 
parameter hypothesis of the turbulent boundary-layer 
literature. Both express the hope that the turbulent 
mixing process can somehow be represented by a single 
empirical relationship describing the response of the 
boundary layer to its environment. Numerous rela- 
tionships, sometimes supported by physical arguments 
and sometimes not, have been proposed to serve this 
need in various engineering applications. However, 
the law of the wake and especially the concept of 
equilibrium flow are new elements in the problem which 
suggest that the office of fourth equation cannot use 
fully be filled until these elements have been accounted 


for. 
(2) KINEMATIC SIMILARIT \ 


a) The Continuity Equation 


In spite of the difficulties mentioned in the Introduc- 
tion, there is one important area of study in which the 
inverse problem for laminar flow can serve as a model 
for turbulent flow. This area involves the idea of 
equilibrium flow, an idea which is not confined to the 
turbulent boundary layer if the term “equilibrium” is 
interpreted to mean invariance of the similarity laws 
for the mean-velocity profile. A restricted class of 


two-dimensional laminar boundary-layer flows with 
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this property was first studied by Falkner and Skan,° 
and later investigations have included certain effects of 
compressibility, of axial symmetry, and of suction or 
blowing at the wall. Briefly stated, the purpose of the 
present paper is to point out and to exploit an analogy 
between the laminar equilibrium flows of Falkner and 
Skan and the turbulent equilibrium flows of Clauser. 
The substance of the analogy will be found to consist 
largely of consequences, inferences, or conjectures 
drawn from the continuity equation and from the in 
formation about equilibrium turbulent flows which is 
available in the law of the wall and the momentum- 
defect law. 

The first of these consequences has already been 
recorded celsewhere;’ it is that the continuity equation 
Ou Ox + Ov OY Q and the law of the wall uw uw, 
f(yu, v) together require the mean streamline slope 
v uw for fixed x to increase linearly with y in some 
neighborhood of the wall. 

A second consequence of the continuity equation 
follows on making a weak and universally accepted 
approximation of boundary-layer type for the external 
flow. If the velocity uw is taken as sensibly inde- 
pendent of y outside the boundary layer and denoted 
by (x), then Ov Oy in this region is a function only of 
v, and the streamline slope v u = v m for fixed x is 
once again a linear function of y. 

These observations mean physically that the diver- 
gence or convergence of the streamlines outside a 
boundary layer is determined entirely by the distribu- 
tion of surface pressure. At the same time, the 
divergence or convergence of the mean streamlines 
within the boundary layer, at least sufficiently close to 
the wall, is determined entirely by the distribution of 
surface shearing stress. These results will be combined 
for the sake of both generality and brevity in the 
form of a theorem, which follows almost automatically 
on introducing appropriate nondimensional variables. 


THEOREM: Consider a shear flow with mean- 
velocity components u(x, y, ¢) and v(x, vy, ¢) such that 
Ou Ox + Ov OV = VU. Assume wu = v = Oat y = Oand 
lim Ou Oy = 7,(x, t) Where uw is a constant. Suppose 
v0 


that uw is independent of y for y larger than some value 


6(x,t). Define u,(x,/) = V7, p where pisa constant; 
. 


u(x, t) = u(x, 6, t); 6*(x, t) = fod — u/um)dy; 1/X = 


—(1/u,)Ou,/Ox; v = wp; and 
O In m, Ov : 
Oln u, Ox , 
6* O ln (6*u,/v) Ox 


Pie. f) =: — (6) 
6 Olnm (m/u,) Ox 


Then the curve obtained by plotting Av Péu against 
y Pé for fixed x and ¢ must leave the origin with slope 
unity and must coincide for y>6 with the straight line 
passing with slope D through the point (1, 1). 


To prove the theorem stated above, note first that 
the formula uu, = f(yu, v) for the velocity profile is 
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equivalent® to the relationship wv u vr. Th 
assumption of Newtonian friction at the wall, becaug 
it requires “ Wu vu, v + higher order terms 


therefore sufficient to establish the result of the theorey 
at the surface. 
Outside the shear flow, on the other hand, u 


u(x,t), and Ov Oy Om;(x, tf) Ox, so that o(a 
v(x, 6, ¢ —(y — 6)Ou, Ox, and, therefore 
Av 6u NV 61) + (Ut, My) |(OM, Ox) (Ou, Ox)| X 
’ 6 - 

where 27, = 7(x, / w(x, 6, ¢). Now in the coordinat 
system (Av 6u, y 6), the straight line defined by the last 
equation intersects the straight line Av 6u y 6 at the 
point (P, P), where 

Fr Av, Oty D\| (1 D 8 
and D is defined by Eq. (5). Eliminating dz éu 


between Eqs. (7) and (S), the integral of the continuity 


equation outside the shear flow becomes 
(1 P)(\v 6u) — 1 = DI’ P)(y 6) — 1 


Finally, from the usual definition of displacement 


thickness 6* it can be shown that 


1 2;)(Ou,;6* Ov (Vy Wy) + (6 y)(OU, On 


Substitution for 7; “; in Eq. (S) yields Eq. (6) for P 
and the result stated in the theorem follows. Not 
that the flow has not been specified to be either laminar 
or turbulent, the mean continuity equation in tw 
dimensions being the same for both, and that the mear 
flow has not been required to be independent of time 


b) Laminar Equilibrium Flow 


The two parameters D and P have so far been de 
fined by completely kinematic considerations. [| now 
propose to show that these two parameters are sepa 
rately constant—1.e., independent of x»—in any one 
of the laminar flows studied by Falkner and Skan, and 
that at least the parameter ? is independent of x in ar 
equilibrium turbulent flow as defined by Clauser. 

Consider first the case of laminar flow, for which the 
statement of equilibrium—i.e., of similarity in_ the 


velocity profile—takes the form 


uu = ¢(y 6) or (mM — u) m = Ply 46 d 


Noting that the shearing stress at the wall is 7 
(uOu OV), = pu,” by definition, it follows from Eq 


) 


that 6u,* «, = constant. Furthermore, the displace- 
ment thickness 6* is evidently proportional to th 
boundary-layer thickness 6. Thus, 

6*, v)(u."/u,") = constant 10 


Substituting for 6*1; v in the definition (6) for P, it 1s 
found immediately that 

6* din (6*i v) 6* » 

P= F = 2? = constant it 

6 dln (um 4, 6 

Turning now to the parameter D, note first that the 

momentum equation evaluated at the wall implies 
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Fic. 1. Kinematie similarity for the laminar Falkner-Skan flows 
u(0°u Oy"), = dp dx = —piuidu, dx. Using Eq. (9 
for u(y), then 
(6°/ v)(du,/dx) = constant (12 


The two conditions (10) and (12) may be combined 


with the two specifications for the variation of 4; with x 


which are found by more elaborate analyses to be 


appropriate for the Falkner-Skan flows. Suppose first 


that 
“uy = Ax" 
Then 6 = Kox' 7", wu, = Kax°"~?4, and 
din uy tn ; 
= = - = constant (13 
dinu, 3n—1 
Alternatively, suppose that 
MW = Cy 
Then 6 Coe~ 8? yn = Cye?™*, and 
din | 
D = = - = constant (14) 
din u, 3 


The exponential velocity variation thus replaces the 
power-law variation when m = +, 

Having identified the parameters D and ? and the 
associated coordinates Av Péu and y Pé, it can be 
Shown that each of the equilibrium laminar flows of 
Falkner and Skan is represented in these coordinates 
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by a single curve independent of Reynolds Number 
Several typical curves, computed from Smith's tables 
for various values of the parameter mw in m4, 

Noting that Pé 26* 


for these flows, it is clear that the curves in the Figure 


constant, are shown in Fig. 1. 


do not depend on the definition adopted for the bound 
ary-layer thickness 6. 

Finally, the fact that the parameters D and P are 
constant for the Falkner-Skan flows implies a relation 
ship D(P) between these two parameters. This rela 
tionship, however, does depend on the definition of 6 
If 6 is arbitrarily taken as the value of y for which 
uu; = 0.99, the curve plotted in Fig. 2 is obtained 
from the numerical solutions tabulated by Smith. 


c¢) Turbulent Equilibrium Flow 


The results just derived for laminar flow are some 
what academic, as a great variety of constant param 
eters can be found to identify the various Falkner-Skan 
solutions. However, consider now the case of equi 
librium turbulent flow, for which the mean-velocity 
profile obeys the defect law 


(uw; — u)/u, = F(yv/6 1) 


outside the sublayer, whereas the local friction law has 
the form 


K(%1/u,) = In (6u,/v) + constant (16 


A calculation of the displacement thickness 6*, neglect 
ing the failure of the defect law in the sublayer,? gives 


*1 
6*u, bu, = | F(y 6)d(yv 6) = constant (17 
ef 
Consequently, 


7 If the true profil given by the law of the wall had been used 
instead of the logarithmic approximation in the sublayer, there 


would have been obtained instead of Eq. (17) the relationship 


6*u,/6u, = constant + |e/(6u,/1 


where the constant € is small (of the order of 75) compared to 


6*u,/v for flow at large Reynolds Numbers. The conclusion that 
the parameter P is constant and equal to 1 7 for turbulent 
equilibrium flow is unchanged when this more accurate expression 


is substituted for Fq. (17) of the text 
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Fic. 2. The theoretical function D( P) for the laminar Falkner 


Skan flows 
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Fic 3. Test of the hypothesis D = d In u,/d In u; = constant 
for equilibrium turbulent flow 
din (6*1 v) u; d In (6u, v) uy 
d In (i u,) 7 u, (wy U,) sia Uu, 
and 
6* d In (6*1 v) 6*u; 
P=-— = Kk — = constant (1S) 
6 dln (mm u,) Ou, 


This expression should be compared to the correspond- 
ing Eq. (11) for laminar flow. For the particular 
defect law given by Eq. (3), moreover, the constant in 


Eq. (18) can be written explicitly as 
P=1+ fr (19) 


This last result, first obtained in October, 1954, 
provides almost the entire motivation for the present 
analysis. For note that the parameters D and P were 
originally defined and connected by one of the most 
general kinematic relationships which might be de- 
veloped for flows of boundary-layer type without being 
so general as to be useless. The relationships (15) and 
(16), on the other hand, are of the essence of contempo- 
rary empirical knowledge of phenomena in turbulent 
boundary layers. The emphasis is on the word 
‘empirical.’ That these apparently unrelated lines of 
investigation are found to converge in the simple Eq. 
(19) must be either a remarkable coincidence or a 
spontaneous manifestation of a fundamental order in 
the problem being studied here. 


d) Kinematic Similarity 


I have shown that the two parameters D and P 
occur naturally together in the integrated continuity 
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TABLE | 


Summary of Experimental Data on D and FP in Turbulent 
Equilibrium Flow 


Reference or = d In u 
Remarks rit? os d In ur D (sevised 
Pure wall flow l ] 
Ludwieg and Till 
mann,” Channel 
VII 1.20 0.02 1.22 +0.03 
Bauer,® 20° slope 1.22+0.02 1.44220.03 
40° slope 1.24+0.02 1.3320.03 
60° slope 1.23 +0.02 1.30 +0.08 
Wieghardt,’ con 
stant pressure 1.55 +0.01 0 
Clauser,* Series 1 92.542+0.05 0O0.862+0.05 0.795 +' 
Series 2 $9320.10 0.752+0.05 O.863 +0 9 


Pure wake flow l 


equation; that D and P are separately constant for any 
one of the laminar Falkner-Skan flows; and that P js 
constant for an equilibrium turbulent flow. These 
observations inevitably suggest the particular hypoth 
esis 


D = (d In m)/(d In u,) = constant 


for the equilibrium turbulent boundary layer. From 
the definition it is seen that the parameter D depends 
on the relative magnitude and rate of change of 1 and 
u,, or alternatively of the free-stream dynamic pressure 
gq and the surface shearing stress 7,. The situation is 
therefore a happy one, in that the hypothesis D = 
constant can be tested a priori by an examination oj 
experimental data. 

Six equilibrium or near-equilibrium flows, omitting 
the case “, = constant, are included in an extensive 
survey of the experimental literature reported else 
where.‘ For these six flows, Fig. 3 shows a plot in 
logarithmic coordinates of m, ™, against uw, u,, where 
u, is inferred from the law of the wall and wm, and u, 
are arbitrary constant reference velocities. The sup- 
position that the data might define a straight line is 
borne out in each case, and the corresponding values of 
D and P are listed in Table 1 and plotted in Fig. 4. 

Hereafter, the term kinematic similarity will be used 
to refer to the property which is implied for an equi- 
librium flow by the statement D = constant. The 
analogy mentioned in the Introduction is seen to emerge 
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Fic. 4. The experimental function D(P) for equilibrium turbu 
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ON THE EQUILIBRIUM 
from this property, inasmuch as the turbulent equi- 
librium flows of Clauser, like the laminar equilibrium 
flows of Falkner and Skan, are then characterized by a 
iw relationship between the two parameters 4; 

This laminar-turbulent analogy will be dis- 
following a 


power-! 
and u 
cussed at some length in Section (4 
demonstration that the 


similarity leads to a solution of the inverse problem for 


assumption of kinematic 


the equilibrium turbulent boundary layer. 


3) SOLUTION OF THE INVERSE PROBLEM 


a) The Karman Momentum-Integral Equation 


Clauser’s main purpose in reference 3 was to show 
experimentally the existence and properties of equi- 


flows with non- 


turbulent boundary-layer 


librium 
vanishing pressure gradient. The appropriate pressure 
distribution was determined in each case by repeated 
trials, working in the direction of the flow. I will now 
show that the assumption D = constant allows Clau- 
ser’s ad hoc determination of pressure distribution to be 
replaced, at least in principle, by an analytical calcula- 
tion 

Three of the four equations needed in the formulation 
of the inverse problem for equilibrium turbulent flow 


are the stipulation 
7 = constant (20) 
together with the local friction law (2 
k(u,/u,) = In (6u,/v) + xc + 22 
ind the momentum-integral equation of Karman, 


Ty p = Uz” = (d/dx)u°6 + 16*(du,/dx) (21) 

For the present, only these three equations will be 
considered, with the four quantities ™, u,, 6, and 7 
taken as dependent variables and x as independent 
variable. 
express the displacement thickness 6* and the momen- 
tum thickness @ in Eq. (21) in terms of the primary 
For this purpose, 


As a preliminary operation it is necessary to 


variables of the problem. as re- 
marked earlier, the concept of a defect law is sufficient. 
In what follows, however, I will retain the notation of 
Eq. (1) in order to show explicitly the dependence of 
various quantities on the single parameter—e.g., 7 
which characterizes an equilibrium flow. 


Substituting the defect law given by Eq. (3) in the 
definitions 
oF = } [1 — (w wm) ldy 22) 
7 
La) 
j= } (u/u)[1 — (u/m) Jdy (23 
70 


and neglecting the small error in the sublayer, the two 
relationships are obtained. 
T Clauser® uses the notation A/é for the quantity called here 


GA/é for the quantity 


r, TE spectiv ely 


‘h(r)/x, and the notation 220( a) /K? 


Rotta™ uses 6*u,/6u- and J,6*u,/6u 


TURBULENT 


BOUNDARY LAYER 5U1 
- 
O° Uu - % 
c= l+f7 ~ dw = (9 4 
ou, J O 
A 6* — @\ u ey \ 
2% =Q+7 In= dw + 
r- 0 i /70 0 0 


In applying these relationships, a convenient first 


and (24 


FD he 


step is the use of Eqs. (25 
Eq. (21) 
Remembering that 7 is constant, the result is 


successively in 


to eliminate @ and 6* in favor of m4, u,, and 6 


A second step is the suppression of one of the three 
derivatives in the last equation with the aid of the local 
form. The natural 


friction law (2) in differentiated 


choice for elimination is d6 dx, yielding 


: 6 duy 1] Hy) 
K kK QQ, — 2(k (Qo — QQ) - 
uy, dx | Hu, Uu, a 


6 du, it) 4) 
K QQ, — 2[k Q. + 20 26 
u. dx u u, a 
A third step, suggested by Eq. (26) as well as by 
experience with the special case “, = constant, is the 


recognition of the quantity 


S = K\M, U,) 


as a fundamental variable. Differentiating the last 


expression, then 


(1/s)(dz/dx) = [(D — 1) D](1/m)(diu, dx) = 
(D — 1)(1 u,)(du, dx) 
where D = d In md In u, by definition. Eq. (26) 


becomes finally 


dx dz D ) 
i — = 2°, — 22, — 2 + 


»)| i 


To recapitulate, Eq. (27) is the momentum-integral 
Eq. (21) evaluated for the special mean-velocity profile 
(1) with constant 7. The defect law (3) has been taken 
to apply throughout the flow, including the sublayer. 
Because no assumption has yet been made concerning 
the parameter D, the effect of the manipulation just 
carried out has been to change the nature but not the 
number of the variables in the problem. It 
conceivable that the form of Eq. (27) would eventually 
suggest the assumption D = constant as a heuristic 


is even 


measure, if attention had not been attracted to this 
hypothesis by consideration of the continuity equation. 
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Note that the flow with D = 0 or m, = constant appears 
as a special case of Eq. (27) for which an integral may 


be obtained immediately, using Eq. (2) to eliminate 6 
in favor of zs. The integral in question was first 
obtained by Karman and is the one cited in the intro 
duction as the only previously known solution of the 
inverse problem for the turbulent boundary layer. 


b) Integration of the Karman Equation 


[ will now assume that an equilibrium turbulent 
boundary layer has the property of kinematic simi- 
larity; 

D = (u./u)(du,/du,) = constant (28) 


Eq. (28) evidently provides the fourth relationship 
needed in the formulation of the inverse problem, and 
is to be considered jointly with the momentum-integral 
Eq. (27), the local friction law (2), and the equilibrium 
condition (20). Now Eq. (28) is itself a differential 
equation which may be integrated immediately, with 
the result 


i / th, = (u./U.,)” (29) 


where m4, = (Xo), U,, = U,(Xo), and Xo is any con- 


venient reference point. Returning to Eq. (27), the 
factor 6 on the left-hand side can be eliminated by 


observing, in view of Eqs. (2) and (29), that 


e /e 1/(1—D) 2-20 
0 09 = (23 Zo) e 
where 6) = 6(x)) and g = 2(x»). Substituting for 6 in 


Eq. (27) there is obtained finally 


~ 
2 


(1/(i-D 1 dg 220, — 220. — Is + 


a 
20% (30) 
4 


1—D 


This equation may be integrated in closed form if 
1 (1 — D) is an integer or half-integer. For example, 


if D = 4/3, 


“" - = (5Q; — 2Q)ao%e *[Ei(s) — Ei(z)] — 
6 
6 ba 
(5Q) — 202) Zo" 2 — | om 
00 20 
6 Zz 6 
(4Q, — 22) ) Zo a —. l cs ag 20. m —_ l ( s1) 
00 So Oy 
Zz e 
where E1(z) = dt 
o t 
and 6/690 = (2/20)~*e" ~*~ (32) 


The quantities (x, 6, z) and (xp, 69, 39) in these expres- 
sions may be considered as variables and parameters, 
respectively, or vice versa, depending on the applica- 
tion. Note that two independent constants of integra- 
tion. 6) and , are encountered in integrating the system 
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Because z = xm u, determines 6y_ ; 


(27) and (28). 


for an equilibrium flow by virtue of the local frictio 


law (2), this means that the two physical scales 6 an 
v/u, may be specified independently at any one statior 

I have evaluated Eqs. (31) and (32) for D f 
numerically, taking = 0.24. These values are in 
tended to represent the spillway flow with 40° slop, 
studied by Bauer® (see Fig. S of reference 4). The 


functions 2;(7) and Q:(7) have been given the forn 
Q = 1+ wand ® = 1 + 1.6007 + 0.7617? require 
by the particular wake function w(y 6) tabulated i 
reference +; the experimental values are not signifi 
cantly different. Similar calculations have been mad 
for r = 0.55 and D = O, corresponding to the flow wit! 
a constant external velocity of 33.0 meters per se 
studied by Wieghardt® (see Fig. 4+ of reference +), and 
for r = 1.54 and D = 4 5, corresponding to the flow 
with moderately rising pressure? studied by Clauser 
(see Fig. 15 of reference 4 Taking the constants § 
and z in each case from the experimental data at a 
point well downstream, the calculated and measured 
values for 6*(x), w:(¥), and 3(v) are compared in Fig 


c) Uniqueness 


One important consequence of the hypothesis D = 
constant for an equilibrium turbulent boundary layer 
is an explicit statement of uniqueness for such flows 
Given a constant value for D, consider first an integral 
of Eq. (30) in the form (x — x0) ‘69 = E(z, 2, a, D 
On eliminating the thickness 6) in favor of 6, 2, and z 
this integral can be written (x — x9) 6 = F(z, 20, 7, D 
Next, note that if P (or 7) and D are separately con 
stant for an equilibrium flow, these quantities must bi 
related by some function D(P) like the one for laminar 
flow shown in Fig. 2. This being so, the dependence oi 
the flow on the parameter D need not be stated ex 


7 The value D = 4/5 (actually 0.795) attributed here 
Clauser’s flow is different from the experimental value D = (0.8 
listed in Table 1. The need for some revision in the origina 
value of D can be argued from the momentum-integral equatior 


in the differentiated form (26 Let this equation be rewritten 


23([0.D + (1 D)| + 2%, 


K°(A/6) = 372,(1 D 


Now the quantity «2/6 can also be evaluated by differentiating 
the local friction law (2) to obtain 


K2(\/6) = [x?/(d6/dx)|[s(1 D) +1 


Elimination of \/6 allows any one of the three quantities 3, / 
and dé/dx to be computed if + and the other two quantities at 
known 
tum balance as well as similarity in the mean-velocity profile 1 


The first equation assumes two-dimensional momet 


the sense required by the law of the wall and the defect law 


The second equation assumes similarity only. Taking + = 1.4 
and D = ().86 from Table 1 for Clauser’s first series of exper 
ments, together with a typical value s = 13, it is found that tl 


computed value of dé/dx is negative. I have therefore reverse 
the calculation, estimating dé /dx from the experiments and calcu 
lating D instead The revised values for D are listed in Table | 
together with an estimate of probable error corresponding to a! 
uncertainty of +20 per cent in dé/dx. That the revisions shoul 
be in opposite directions for Clauser’s two flows is suggested ! 
slight discrepancies in momentum balance noted elsewhere (¢! 


the functions #(x) in Figs. 15 and 16 of reference 4 
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plicitly; thus (v — 4%) 6 = G(s, %, @ Finally, 
suppose that a convention can be found for specifying 
the origin or initial point x», 6), 2, ete. in such a way that 
depends on am alone. For example, require the 
momentum thickness # to vanish for 4 \ 0, and 
note that Eqs. (24) and (25) then imply «am u Z 
20.( ar) Qa ar rhe integral of Eq. (30) under 
these conditions can thus be expressed ultimately in the 
irl (2, 3 But if 6u, v and 6 x for constant 


rare tunctions Ol KH; U. alone, then so are ux v 


} So also are 6” 6* 6. 4,6 v, and similar 


uantities, by virtue of various relationships already 


lerived from the momentum-defect law 
therefore, the 


Given the property D constant 


onclusion is that quantities like 6* 6, Cy 2u,?/ Uy 
ind J< u,9 v can be expressed as one-parameter 
functions of a uniquely defined streamwise Reynolds 
Number wx v for the class of equilibrium turbu- 
lent flows even when wu; depends on x. This conclusion 
does not require the assumption that the wake function 
y 6) in Eq. (1) is universal, because the mean 
velocity profile in an equilibrium flow is adequately 
f this argument by the 


expressed for the purpose « 
defect law (4 Neither would the conclusion stated 
here be changed if the exact mean-velocity profile in 
the sublayer had been considered in Eqs. (24) and (25), 
is the quantities Q)(7) and Q2.(7) could then be replaced 


by Q:(7, 2) and Qe( 7, 2 


d\ The Shearing-Stress Profile 


The distribution of shearing stress within an equi 
librium turbulent boundary layer may be found by 


integrating the boundary-layer equations 


(Ow Ox) + (Ov OV) = O 


u(Ou ON v(Ou Oy uy(du, dx) + 


1 p)(Or Oy 34) 


for the profile of Eq. (1) with constant 7. This caleu 
lation, which amounts to inferring the forces applied to 
the fluid from the observed accelerations, should be 
equivalent within the boundary-layer approximation 
Noting 


by definition, the result of a tedious lot 


to a direct measurement of the shearing stress. 
that 7 pl 


of algebra is 


\ 6) 1 tT Te) |] = (y 6) [€? — Dz? + (1 D)ts X 
l 21 — D S\ Wyo — WwW -_ 2F a T 2w Sh 
where 2 KM, U,,¢ = ku U,, and a, and a» are incom- 


plete integrals corresponding to 2; and ®: in Eqs. (24 
and (25 The functions w)(7, vy 6) and we(7, y 6) are 
defined and tabulated in reference 4. 

5,000 for four equilibrium flows 
which have been observed experimentally—first, P = 
24, D [.33: second, P = 1.55, D = 0; third, 
P = 2.54, D = 0.80; and fourth, P = 4.93, D = 0.86 


t 


the mean-velocity profile according to Eq. (1) and the 


Taking 6u, v 


total shearing-stress profile according to Eq. (35) are 
shown in Fig. 6. Also shown as a crosshatched region 


is the velocity defect in the equivalent wake.‘ 
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either class of flows can be represented by a certain 
function D(P). To study this analogy more closely, 
the flows in question may be arranged according to 
increasing values of the parameter /, as in Figs. 2 and 4. 
Noting that the present analysis of turbulent flow is 
based on a defect law with four characteristic param- 
eters ™, “,, 6, and P, two points of coincidence appear 
at D = Oand D = o, 
constant and u, = constant, respectively. 
6 = constant, which in the laminar Falkner-Skan 
family is the plane stagnation-point flow, apparently 
does not appear among the turbulent equilibrium flows. 


corresponding to the cases 4; = 
The case 


Two other special cases will be investigated here as 
possible points of coincidence in the laminar-turbulent 
The first is the case uv, = 0, and the second 
Both will be treated by 


analogy. 
is the case “4, u, = constant. 
slight extrapolations of the ideas in Sections (1) and 
(2), and the arguments given here may be accepted or 
rejected on their merits without seriously prejudicing 
any of the earlier discussion. 

The profile parameter 7 in Eqs. (1) and (2) is a 
measure of the relative magnitude of the wake and wall 
components in the mean-velocity profile. According to 
Eq. (2), m must become indefinitely large if u, ap- 
proaches zero while 6 and m remain finite. At the 
same time the product 27u, xu; approaches unity, so 


that the mean-velocity profile for 7, = 0 is locally 
u/u, = (1/2)w(y 6 (36) 
where w = 2 when y = 6 by definition. 


I will now define a pure wake flow as an equilibrium 
turbulent boundary layer for which the surface shearing 
stress is zero everywhere. The profile is then given by 
Eq. (36), the parameter P = 
quantities 6* 6, 6/6, and 6* @ are constants. 
be expected, the laminar equilibrium boundary layer 
with 7, = Ois one of the Falkner-Skan family—namely, 
the flow for which P has its largest value in Fig. 2. The 
velocity profile for this laminar flow can also be repre- 
sented by Eq. (36), although with a different function 


1 + 7 is infinite, and the 
As might 


w(y 6). 
Consider therefore the momentum-integral equation 
of Karman, Eq. (21), when 6 6 = constant and 7, = 0; 


( 4 ( .\ m 
2+ = QJ (37) 
6 iy dx 


+ 6*/0 


ee 
Ou, = 


1 dé 


6 dx 


thus, 


It follows that constant 


The shearing-stress profile is readily obtained by 
integrating the boundary-layer Eqs. (33) and (34) 
directly for the velocity distribution (36). The result 


can be written 


vu y es *) fe 
: == 2 wd¢ 
dé dx 6 2+ 6* 6/ wo 
( l )? C5 yep . e+ 
= = war 
2+ 6*/o/6 \2+6*/0/4J0 
6* 6 I y/b 
a: ae { wedge 
2+ 08 (OF ta 


(39) 


T/ puy- 


dé dx 


(40) 


1 


1OoSs7 


ICAL SCIENCES TULY, 
If the flow is laminar, the shearing stress 7 in Eg 
(40) must be related to the velocity by the Newtonia; 


wou Oy = Substituting jy 
and noting that the left side, like the right 


law 7 = (ut, /26)w’(y/d). 
Eq. (40 
must be a function of y 6 only, it follows for the con 


tinuously separating laminar equilibrium flow that 


(6u, v)(d6 dx) = constant 


Combining this condition with Eq. (38), it is seen that 


nx /p 


1/(3+26 
UX = constant 4] 


where it has been assumed that “; = « 
expression has the proper form m.x~" = constant, an 
the value of the exponent is known’ to be very nearh 
n = —0.0904 for the flow in question. In the present 
context, however, it is instructive to form an inde 
pendent estimate of this exponent by taking 6* @ = 
a value consistent with the expected antisymmetry oj 
the velocity profile. Then » = —1/(3 + 26*/6) = 
—1/11 = —0.0909. 

Next consider turbulent flow, for which the Newto- 
nian stress law r = wOu Ov is no longer valid and has t 
be replaced by some other assumption. Taking int 
account the condition 7, = 0, it is reasonable to requir 
the mean motion to be independent of the viscosity, as 
in the case of free shear flow. If so, then only tw 
lengths 6 and x can appear in the problem, and thes 
two lengths must be proportional. But if dé dx isa 
constant, it follows from Eqs. (39) and (40) that 1 
and 7 q are functions of y 6 only; and the converse is 
also true. Finally, replacing 6 in Eq. (38) by x, there is 
obtained 


1/(2+6* 
UX = constant (42 


If 6* 6 is again taken as 4, the exponent is 1 ‘6, whik 
the supposedly more accurate value 6* @ = 4.18 from 
reference + gives 0.162. In any event, the difference 
between the exponents in Eqs. (41) and (42) for laminar 
and turbulent flow is consistent with the empirical 
observation that turbulent flow will in general support a 
more rapid pressure rise. 

Finally, it can be shown that the parameter D must 
be equal to unity for the turbulent pure wake flow 
Eliminating dé dx between Eqs. (37) and (39) with the 
6, substituting for 7 a in Eq 
—(X/u,)du; dx, then 


latter evaluated at y = 
(S) for P, and noting that D = 


(2)4) 
Hi -~ Dp 6 ‘ 


If the flow is laminar, then P = 26* 


unity and 6* 6 must be nearly 4, so that (1 — D 


P = 


6 must be near} 


De-= 


5/2and D = 27. The exponent » obtained from Eq 
(13) of Section (2) is again —1/11. If the flow is 
turbulent, on the other hand, the parameter P = 1 +7 
is infinite in Eq. (43), and therefore D = 1. 


The mean-velocity distribution (36) and the shearing 
stress distribution (40) are plotted in Fig. 6 for the 
particular wake function w(y/6) of reference 4. The 
shearing stress being computed as 7 Taz, it is not 


atx =0. This 
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ON THE EQUILIBRIUM T 


] 


necessary to specify a value of dé dx for the flow in 


questio! 


b) The Pure Wall Flow 


A second special case of considerable interest is pure 


wall flow, obtained from Eq. (1) by putting x = 0—.e., 
The mean-velocity profile is then given by the 


from which it follows® 


P= l. 
law of the wall, wu, = 
that wu, and yu, v are constant on mean streamlines 


f(yu,/v), 


and that 
Av 6u = ¥ 6 + 
The shearing- 


where 1 A = —(1 u,)du, dx as before. 


stress profile is most readily obtained by putting 


= w = 1 in Eq. (35 


6, therefore, where u = m4, and 7 = O 


r uu, \° Uy 
“* = (. ) (1—D) - 2. ) +? 4) 
6 Uu, u, 


A number of arguments can be found, all of them 


At y 


unfortunately somewhat porous, for supposing that if 
the parameter D is constant for the pure wall flow it 
ought to have the value unity. Certainly the point 
P = D = 1 is favorably located in Fig. 4 with respect to 
the experimentally determined points in accelerating 
flow, if the hypothetical function D(P) for equilibrium 
turbulent flow is to resemble the one in Fig. 2 for laminar 
flow. The statement D = 1 can also be argued on 
dimensional grounds; the first of the two characteristic 
lengths 6 and v/u, does not enter the problem explicitly, 
and the two lengths in question should therefore be 
proportional to each other. Finally, the statement 
D = 1 follows if v/u and 7/7, in Eqs. (44) and (45) are 
required to be functions of y 6 alone, and conversely. 

If D = 1 for the pure wall flow, then 6u, v, ™ u,, 
and 6 \ are constants, and, therefore, wx v = constant. 
The external flow is either a source or a sink flow, and 
only the latter is permitted by the conditions \ 6<0 and 

u, < 0 which are obtained from Eqs. (46) and (44) 
The 
corresponding laminar sink flow is readily identified as 
the one for which P takes on its smallest value in Fig. 2. 
This flow may be studied (1) directly;!! (2) by making a 
boundary-layer approximation in the known exact 
solution of the Navier-Stokes equations for radial sink 
flow; or (3) by a suitable limiting process in the Falk- 


When “; and “, are taken as positive quantities. 


ner-Skan equation, letting the exponent ” in m4x7" = 
constant approach — 1 from below so that the Falkner- 
Skan parameter 8 = 2n/(n + 1) approaches +@. 
Because vu = y/x for this laminar flow by assumption, 
the friction coefficient must be independent of x and the 
velocity profile must be expressible® in the form uu, = 


g(yu,/v) with u,/u, aS parameter. For a given 
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laminar flow there is only one function ¢g(yu, v) such 
that the boundary-layer equations are satisfied at the 
same time that r = wou Oy. For turbulent flow, on 
the other hand, the function ¢(yu, v 
in the law of the wall, 


has to be identi 
fied with the function /(yu, v 
and the momentum equation is automatically satisfied 
is computed there 


because u(x, Vv) is given and r(x, y 


from. In either case the friction coefficient is constant 
but its value can apparently be chosen arbitrarily. 
The 
corresponding shearing-stress profile in the pure wall 
= 5,000 and D = 1. 


for values of 


logarithmic mean-velocity profile and the 
flow are shown in Fig. 6 for 6u, v 
I have also computed the function 7 7 
6u, v ranging from 10° to 10°, and have found that the 
various curves can scarcely be distinguished in the 
Figure. That is, r 7, as a function of y 6 is for prac 
tical purposes independent of the Reynolds Number 
6u, v for the pure wall flow. These calculations refer 
to physically different boundary layers, not to different 
stations in the same boundary layer. 

Finally, if the conjecture D = 1 is in fact correct for 
the pure wall flow, then the edge of the boundary layer 
is a mean streamline and the net velocity of propagation 
of the laminar-turbulent interface into the nonturbulent 


fluid is zero. 


c) An Interpolation Scheme 


In any practical application of the concept of kine- 
matic similarity—for example, to diffuser design—some 
interpolation method is needed to supplement the 
The 
method proposed here is far from unique, and depends 
on the development of two quantities which, unlike the 


experimental values of D(P) listed in Table 1. 


parameters D and P, remain finite as 7 increases from 
One such quantity is the strength of 
zs. Another is 


zero to infinity. 
the equivalent wake, 27u, xu, = 27 
the rate of mass entrainment in the boundary layer. 


Defining 


(dé dx) — (0) / 8s = 
(1 pu;)(d/dx)[puy;(6 — 6*)] (47) 


(§ uy) = 


then the quantity s is seen to be the net velocity of 
propagation of the boundary y = 6 with respect to the 
free stream. Now the ratio su, may be expressed in 
terms of the local friction coefficient and the parameter 
7 with the aid of Eqs. (2) and (18): 


su, = (6 s)(dz dx)(s — mr) 


where z = xm, u,. But (6 s)ds dx is a known function 


of z, 7, and D from Eq. (27); thus, 
g7)(1 — D) 


ma (48) 


Je) 


The case of pure wake flow is most conveniently treated 


separately to obtain 27 z = 1, and, from Eqs. (47) and 


= ] ee) 
«2, = - 2x? dx\2 + 6*/0 


(39), 
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For given values of 7 and D and for a specified value 
S K°u; may be 


Plotting s Ku 


of 6u, v, the quantities 2, 27 and 


computed from Eqs. (2) and (48). 


against 27 z for the six flows of Fig. 6, estimating 
dé dx = 0.252 for the pure wake flow from reference 4, 
and using essentially straight line interpolation except 
near 7 = O, the inverse calculation? for D(P) leads for 
6u,/v = 5,000 to the curve shown in Fig. 4+. This cal 


culation, as might be expected, is not at all sensitive to 


the value chosen for 6u,/ v. 


(d) The Hypothetical Function D(P) 


Perhaps the most instructive physical interpretation 
of the new parameter D comes from the fact that the 
mean streamlines in a shear flow intersect at a common 
origin for any region in which v w for fixed x is a linear 
function of y. For a boundary layer one such region is 
the one near the wall, including the sublayer in the case 
of turbulent flow; herev wu = y XA. In the absence of 
viscosity there is a corresponding relationship vu = 
Dy for the ambient flow. Thus the parameter D is a 
measure of the extent to which any divergence or con 
vergence of the streamlines in the external flow, which 
is to say any pressure gradient, affects the streamlines 
in the neighborhood of the wall. This interpretation 
applies equally for laminar and turbulent flow under 
quite general conditions, and like the defect law itself 
does not involve the viscosity of the fluid explicitly. 

At the beginning of this paper the inverse problem 
was defined as the determination of the surface shearing- 
stress distribution 7,(x), given the ambient pressure 
distribution p(x) and suitable initial and boundary 
conditions. For a general flow a direct solution other 
than an experimental one is not yet possible for this 
problem. However, a knowledge of p(..) is sufficient to 
determine the streamline pattern outside the boundary 
layer, while a knowledge of 7,.() is sufficient to deter- 
mine the streamline pattern near the wall. For equi- 
librium flow the assumption of kinematic similarity 
connects these two observations directly, thus providing 
a detour as it were around the obstacle represented by 
the turbulent transport mechanism. The connection 
may be viewed either as a relationship between stream- 
line patterns in two regions of the flow or simply as a 
power-law relationship between g and 7,.. 

The second of these two interpretations shows more 
clearly the role which is played by the concept of kine- 
matic similarity in the solution of the inverse problem. 
For an equilibrium turbulent boundary layer it is evi 
dent from Eqs. (2) and (26) that the Karman momen- 
tum-integral Eq. (21) can be reduced to the form dx = 
Gi(m, u,, P)du, + Ge(m, u,, P)du,. The assumption 
D= constant then leads to a 


tT For Clauser’s second flow it appears that the parameter D 
ought to be taken as 0.866 if a smooth curve is to be obtained in 


D «<= 


constant or 4H, = 


9. /e 


the coordinates (27/z, s/x?u;) 


AERONAT 


I 


ial SECISBNCES PUL, £9534 


separation of variables, and the only further informa- 
tion required for a full solution is the empirical function 
1G ag 

What is ultimately needed, of course, is not a physical 
interpretation but a physical principle, from which 
might be deduced not only the existence of a function 
D(P) for equilibrium turbulent flow but the form of this 
function. Granted that Fig. 3 justifies the assumption 
D = constant as an interpolation device for the data of 
the Figure, other reasons must be found for making the 
same assumption for a general equilibrium flow. The 
theorem presented earlier, in which the two parameters 
D and P first occur, amounts at best only to circum 
So does the analogy proposed here 
In the 


stantial evidence. 
for laminar and turbulent equilibrium flows. 
absence of a physical principle, therefore, any discussion 
of kinematic similarity requires an act of faith in that 
neither of the two statements D = constant or P = 
constant can be said to imply the other. The problem 
having been stated in these terms, it follows from ex- 
perience with the special but by no means trivial case 
“u, = constant or D = 0 that a further attempt must be 
made to account directly and specifically for the con- 


cept of a defect law 
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A Theoretical and Experimental Study of 
Airplane Dynamics in Large-Disturbance 
Maneuvers 


DONALD W. RHOADS* ann JOHN M. SCHULER** 
Cornell Aeronautical Laboratory, Inc. 


SUMMARY 


This paper reviews an effort to solve some problems of large- 
listurbance maneuvering flight. The period covered is from 1949 
to the present, and results of both theoretical and flight test in- 
ve stigations are presented 

rhe rolling pull-out, a critical maneuver for vertical tail loads, 
was chosen for study Equations of motion were developed, and 
solutions were compared with responses obtained from flight tests 
The equations, validated by this comparison, were then applied to 
critical present-day airplane configurations, and a general study 
of large-disturbance motion was made including variation of per 
tinent parameters 

This paper represents a first step in the study of nonlinear as 
pects of airplane dynamics and control in large-disturbance 
maneuvers. Results show that inertial and aerodynamic non 
linearities introduce destabilizing tendencies not predicted by 
linear theory Although complex equations are required to pre 
dict motions accurately, the characteristics of the motion may be 


described using simple mathematical techniques 


SYMBOLS AND DEFINITIONS 


Che orientation of the airplane body axes with respect to earth 
xes or inertial space, defined by rotations y, 6, and ¢, is shown in 
Fig. | The orientation of wind axes, stability axes, and body 
ixes and the conventions for linear velocities, angular velocities, 
forces, moments, and control surface deflections are shown in 
Fig. 2. The symbols and definitions used in this paper agree 
substantially with those recommended in reference 13 

> = stability axes 


K % = body axes 


rhe body axes (unless otherwise noted) are oriented so as to 
coincide with the stability axes for straight and level flight in the 


reference flight condition 


,9, ¢ = orientation angles of earth to body axes (see 
Fig. | 

p = angle of side slip see Fig 2 

a = angle of attack (see Fig. 2 

= control surface deflection 

= angle of attack of the x principal axis 

= components of linear velocity 

= components of angular velocity 

, = component of acceleration along 2 axis, g’s 
Y,Z = components of force 
Vv, N= components of moment 

\ = load on the vertical tail 

\ = load on the horizontal tail 

” = mass 


= time 


Presented at the Stability and Control Problems During Large 
Maneuvers Session, Twenty-Fifth Annual Meeting, IAS, New 
York, January 28-31, 1957 


“Research Aeronautical Engineer, Flight Research Depart- 


Associate Research Aeronautical Engineer, Flight Research 
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Cp 
Ck, 
Ck, 
Ck 
Cka 


Cxa 


weight 

moments of mertia 

product of inertia 

air density 

altitude 

Mach Number 

total airplane velocity = V W 

reference forward velocity along the x axis 
dynamic pressure, (1/2)pV? 

wing mean aerodynamic chord 

wing span 

acceleration of gravity, 32.2 ft./sec 

aspect ratio 

sweep angle, sweep of quarter chord line 
¥/(1/2)pV2S, Y/(1/2)pV2S, Z/(1/2)pV2S 
L/(1/2)pV2Sb, M/(1/2)pV2Sc, N/(1/2)pV2S! 
Lift /(1/2)pV2S 
Drag ] 2)pV? 
OCK /O( pb /2V K = x,y, 2, l,m,» 
OC K /O(ge/2) 

OCK/O(rb/2V 

OC K /Oac/2) 

OCK/Oa; 


Derivatives with respect to 8, 6¢, 6,, and 6, are defined similarly 


those with re spect toa 


CL; 2 = QC ,/06 
Cms = OC,,/06 
AC,» = 40C,/O(pb/2V = side force due to rolling 
velocity— dihedral contri 
bution 
C,,/a = [0C,/%rb/21 a = variation of rolling moment 
coefficient due to yawing 
velocity with angle of 
attack—-untwisted wing 
contribution 
AC, = AlOC,/O(rb/2) = (), due to wing twist 
( a = [0C,/0(pb/21 @ = variation of yawing moment 
coefficient due to rolling 
velocity with angle of 
attack—-wing contribu 
tion 
AC,, = A[OC,/O( pb /21 = C,,, due to wing twist 
C,,/a? = [0C,,/(rb/21 a? = variation of vawing moment 
coefficient due to yawing 
velocity with angle of at 
tack 
at, = A[OC,,/0(rb/2) = C,, due to profile drag 
Nba: N(éat) = OC,,/06. = fla) = yawing moment due 
to aileron as a func 
tion of angle of at- 
tack 
ly, /4 = distances from the c.g. to the 


center of pressures of the 
vertical tail and horizontal 


tail, respective ly 
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2. 


Orientation of earth and body axes 


Fic. 1. 


= distance from x body axis to center of pressure of ver- 
tical tail 
a4, = angle between the zero lift line and the x body axis 
in reference condition 
Qs = angle of attack for a given 6, in a steady pull-up, index 
of @ in a rolling pull-out 


Pr = roll rate measured at @ = 90°, index of p in a rolling 
pull-out 
Perit = Critical roll rate 
é.¢.7,%.%. 
| pe OS ea 
B, Cc, F,G, H, J,$ = Mach Number and airplane geometric 
R, L, P,Q. R. constants used in this paper 
S, 7, 0, W, 2X 
Dimensions 


time, sec 

distance, ft 

angle, rad. (unless otherwise noted 
force, lb 


moment, ft.-Ib 


mass, slugs 


Subscripts 


e = earth, component or axis 

1 = wind, component or axis 

s = stability, component or axis 
é = elevator 

a = aileron 

r = rudder 








&s 


w 


Orientation of body and stability axes, linear and angular 
velocities, and control deflections 


Fic. 2 


TICAL SCIENCES-—JULY, 1957 


FRL = fuselage reference line 

] = reference conditions 

iW = wing 

We + F = wing plus fuselage 

H = horizontal tail 

V = vertical tail 

0 = incompre ssible e. ¢.. c 


ai 
= derivatives in which all components refer to st 
bility axes—e es Gas = OCn O( p.b/21 


INTRODUCTION 


a. Wor_D War II, investigations by the 
NACA established that certain vertical tail fail- 
ures were caused by maneuvers combining a roll and 
pull-up, simultaneously. Large sideslip angles, result 
ing from the yawing moments, were responsible. Thus 
a need was established for designing structurally for 
such a maneuver, which was named a rolling pull-out 

Simplified equations developed for these early tail 
failures were found to be inadequate for a later airplane 
Theory did not agree with flight test, so more complete 
equations were developed. Further experiments and ef- 
forts by aircraft manufacturers to design for this 
phenomena led to considerable confusion. A number of 
different equations were developed and proposed for 
use in design. 

In 1949, the Structures Branch of the Aircraft Lab 
oratory, WADC, contracted with Cornell Aeronautical 
Laboratory to perform studies and flight tests as sug 
gested by Joseph Kelly, Jr., of WADC. 
plete a set of equations as possible (limited only by the 


First, as com- 


available data and practical details of computation 
was to be written for the rolling pull-out maneuver; a 
typical problem was to be programed using these equa- 
tions; and simplifications were to be determined which 
caused no significant loss in accuracy. Next, work was 
to start on the measurement of tail loads and airplane 
motions in flight tests of the rolling pull-out maneuver. 
Finally, the equations previously developed were t 
be used to calculate motions and tail loads for the pur- 
pose of comparing theory and experiment. The actual 
control inputs measured in flight were to be used 
Solutions for the initial theoretical study were obtained 
by the Institute for Numerical Analysis of the National 
Bureau of Standards.'* Most of the computations for 
the final stage were performed on analog equipment at 
the Aeronautical Research Laboratory, WADC. This 
work, using an F-SOA airplane as a test vehicle, was 
successfully completed after several years. The re 
sulting five degree of freedom equations found consid- 
erable use. 

However, by this time, many investigators were find- 
ing that certain trends in airplane design resulted in 
rolling pull-outs which were so extreme as to be intoler- 
able both structurally and to the pilot. Therefore, in 
1953, under the administration of Philip Parmley and 
Walter Boccius of the Structures Branch, Aircraft Lab- 
oratory, WADC, the Cornell Aeronautical Laboratory 
was placed under contract to study this maneuver as it 
applied to modern fighter aircraft. The effects of 
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changes in airplane configuration, flight conditions, and 
control motions were to be studied. Though the cal- 
culation of tail loads remained the primary objective, 
emphasis was placed on obtaining a fuller understanding 
of the dynamics, since once the motions are satisfac- 
torilv defined, the loads can be predicted with confidence. 

A digest of the results of this series of investigations is 


presented in this paper. 


INITIAL THEORETICAL STUDIES 


Initial theoretical studies of the rolling pull-out ma- 
neuver proceeded with the thought that the equations 
of motion should represent large-disturbance maneuvers 
as completely as possible. This immediately implies 
that coupling between the lateral and longitudinal 
modes must be included, since, by definition, separation 
based on small-disturbance linearized theory is no 
longer valid. 
coupled equations of 


This concept of “complete” 


motion requires important inclusions—namely: (1 

Inclusion of an expression for the summation of forces 
along the x axis, or ‘‘forward speed change’ equation, in 
the equation system. This has an additional note of 
complexity in that to be consistent, changes in Mach 
Number throughout the equations must be taken into 
disturbances 


account. (2) Inclusion of products of 


(such as wp, vr, and gr). (3) Inclusion of variations in 
the stability derivatives due to angle of attack (also 
Mach Number) and other aerodynamic coupling. 

Because of the practical limitations, several restrict- 
ing assumptions had to be made. 

1) The effects of nonstationary flow are omitted. 
Previous work had shown that this omission will not 
Additional degrees of 
The 


aircraft to which these equations were applied (F-SOA) 


cause significant errors. (2) 


freedom due to aeroelastic effects are neglected. 


has sufficient structural stiffness so that structural vibra- 
tion frequencies are far removed from those of airplane 
motion as a rigid body. The systems may be considered 


independent, and only quasi-static effects need be in- 
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cluded if required. (3) Engine gyroscopic effects were 
felt to be small for the configuration under study and 
were not included. 

The system of equations thus derived represents a set 
of six simultaneous nonlinear equations with variable 
coefficients, describing motion along and about all three 
axes of the aircraft. 

Animportant step toward the practical application of 
the equations of motion to the prediction of motions 
and loads for a new aircraft design was that of simpli- 
fication. It was obvious from the initial considerations 
of the problem that the derived equations would require 
considerable effort to solve, even if use was made of 
available machine techniques. Thus, a plan of simpli- 
fication was devised with the purpose of arriving at a set 
of equations, more adaptable to solution, but without 
undue sacrifice of accuracy. This was considered a 
“broad”’ plan with the thought that the simplified sys- 
tem selected might warrant further study along these 
lines. 

Since the primary unknown effect was that of speed 
change, the obvious procedure was to eliminate this 
equation from the set, and compare the solution with 
that of the complete system. This promised to be a 
rather far-reaching simplification, since Mach Number 
functions, which were included for essentially all the 
derivatives, would be eliminated. A second simplifica 
tion (which also included the first) was to decouple the 
aerodynamic terms. This was thought to be indeed 
drastic, since even cursory examination of the rolling 
pull-out maneuver reveals a dependency of yawing 
moment on angle of attack. A third simplification, 
somewhat academic, was to (including the above two 
simplifications) the This 
latter step would resolve the motion to that of an aileron 
roll, and while not thought to be critical, would provide 


uncouple inertia terms. 


an interesting comparison. 

The complete equations of motion, derived on the 
basis of the assumptions noted at the beginning of this 
section are found below. The inputs are elevator and 
aileron; a rudder locked condition is assumed. 


sak MR + 2 a 
m\i — vr +wq| = —Wsiné+ C., agS sin a — 
RV1—- 1+ 2 
+ i k Cis ( ia. 7+ @ i ia gS cos a + Cp UyS + iu Sin a la 
RV 1 — M242 J 2 
a ” R, 2 S| : 
n sas wp = r(l a ee = 4 13 _ ¢ : Bq — 
rene wo \RrV1 -— 24+ 2/7 9 
Ry +2 pz F 
C ( J ) J \ qSi die J AC, | Ph 
ov \RyV1 — M2/ Ui +4 TGV1— 12+ H) 
l MR + 2 pb ‘ 
ag ( ) (as +a ( - gS + IV cos @ sin @ lb 
RR \RV1 — M2 + 2 Paci +n 
R+ 2 
miw — (lL, + u)g + vp] = W cos @cos ¢ — C; ( 7 ) agS cos @ — 
“\RV1- IP 
R+ 2 ; , —_ . 
yi + fat gS sin a + Cm, gS ; 56. cos a — W cosa le 


) 
E. + ki C, ( : ) (a4, T @ 
*\RV1-— I? 


‘H 








O10 JOURNAL OF THE AERONAUTICAL SCIENCES JULY, 1957 


F \ Ry + 2 by?Sy 
[,p — Inzpq + Uz — 1,)qr — Int = | Coon ( } + C DH ( : ) — 4 x 
GVi1— MW?+ MuV1 — 12+ 2/ 8S 


Ta + Ud, — 1,)pr + I2p? — Ir? = JBC, ( 
te RVI 


+ . Ry + 2 
(Qa4a + a’) ( gsc — ©, 
“Ho \RyV 1 — M2 + 2 


— : U 
(F+ GM + AM) |! 


Te+ Ud, — [,)pq — [22p + Inrq = ge eC i ( 


KR + 


F Ry, a 2 
ie ( ) By Sb _ "ie ( ) BqS\ sy + 
ee NEM I AJ? + H/ = RyV 1 — MWe+ 2 
wevicel 
(a4, +t a) + AC, 
a KV 1 —- IP 
r F 
lyzySvqg —— +C ( ) aug ld 
Ui + u eo \GV1 — M?+H 
2 . R+2 
) e+ c| C, ? ( ) 
M? +2 = RV1 —-— W242 


x 
1,” Ry + 2 
) 4 qSu ‘ine C; ( : ) x 
Ur; + u “H \RuV1 — M?4+ 2 


rb? 54 C ( Ry, si 2 ) x 
2U +0) 7 te a — ae + 2/ * 


T 2 =- 


(U, + u)- 


9) 


WU L?Sy . . : 
r q+ C,, dq le 
d Oe 


ibe ze 


Ry + 2 Syly | ‘ 
: 8qgSb a 
RyV 1 — M2 4+ 2/ Sb 


P pb 


c. . 
| ? (aa, + a) = ——— — =< At x : gS + 
a (RV 1 — 124+ S\(OV1- 124+ T 2.01 + u 


P 
[N(bnp) — NO, 


' eh a oe ae = 
. t (RV 1 — M2 + 5)(0V 1 — M2 + iT) 


In addition to these equations, others were needed to 
describe the motion, and further point to the complexity 
of the system. An equation used to express total nor- 
mal acceleration in gravity units was 


») 


n, = \[w — (Ui + wq + vp)/g} — cos@cos@ (2) 


The following trigonometric functions were used for the 

angle of attack and sideslip: 

a = tan—[w/(Ui + u)] B = tan—! [v/(U; + u)] 
The orientation angles, ¢, 0, and y were expressed as 

= S(p + qsingtan@é+rcos@tan 6@)dt (3a) 

6 =f (q cos @ — rsin ¢)dt (3b) 

¥ =JS[(¢ sin ¢ + rcos ¢)/cos O]dt (3e) 


The tail load equations, which are based simply on the 
motions obtained from the above, were defined as 
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u as used in the equations in this section refers to incre- 





mental forward velocity, and the dynamic pressure 


g = (p/2)((Ui + u)? + v? + w?). 


No attempt will be made here to describe in detail 
the derivation of the various terms in the equations of 
motion; this may be found in reference 1. The accel- 
eration terms (left-hand side) are shown in well known 
form; the aerodynamic terms (right-hand sides), though 
general in overall concept, are based in detail on the 
aerodynamic data available at the time for the F-SOA 
airplane under the initial flight conditions of level flight 
of Mach Number 0.60 and 20,000 ft. altitude. It is 
important to note the aerodynamic nonlinear and 
coupling functions included. The drag component is 
partially a function of sideslip angle (the term /9"); 
the pitching moment is represented as a nonlinear func- 
tion of angle of attack; and the lateral derivatives 
Cup» Cur, and Ons, are functions of a. It is obvious 
that these particular variations pertain to the configura- 
tion under study, and that the nonlinearities must be 
considered separately for each configuration. 


It is obvious from examination of these equations 
that one of the major contributions to their complexity 
is forward speed change, particularly as applied to Mach 
Number variation. Due to complexity, it was neces- 
sary to employ the slower but larger capacity digital 
computer rather than the faster analog type. Solu- 
tions for this initial study were made using the digital 
facilities of the National Bureau of Standards Institute 
for Numerical Analysis at Los Angeles. 
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Comparisons of the solution of the complete equations 
on with those employing the aforementioned 


mplifications are shown in Fig. 3. As analysis of the 





solutio1 


for the complete equations proceeded, it be 
eame obvious that on the basis of the resulting motions, 
additional “local” simplifications could be made. Thus, 
the first simplification (speed change neglected) could 


be extended to include the following 
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It must be reiterated that such simplifications are based 
n the initial conditions and resulting motions; each 


case must be considered separately. 

There are two interesting conclusions which may be 
drawn from this initial study comparison. They are: 
1) For the where derivatives do not 
rapidly with Mach Number, the forward speed change 


case change 
equation (and Mach Number functions) may be neg- 
lected. For flight conditions where the above is not 
so, the effect of speed change should be further inves 
tigated. (2) Uncoupling of the aerodynamic terms re 
sults in unconservative values of sideslip angle and ver 
tical tail loads of such magnitude as to render this 


simplification invalid. 


FLIGHT TEST PROGRAM 


As part of the overall initial rolling pull-out study 
program, a flight test phase was instigated, the pur- 
pose of which was twofold: (1) to provide a wealth of 
actual rolling pull-out flight test data, which although, 
because of practical considerations, remained within 
the safe structural limits of the airplane, could be used 
as a basis for further research on vertical tail loads. 
Through extrapolation or more severe flight tests, the 
critical rolling pull-out maneuvers for vertical tail de- 
sign for structural integrity might be determined. (2 
to use such flight test data in an effort to validate the 
equations of motion. The obvious importance of this 
is that accurate predesign prediction would be possible. 
This plus the determination of the critical rolling pull- 
out maneuvers noted above would greatly aid in the de- 
velopment of vertical tail design criteria. 

An F-SOA airplane was instrumented to provide mo 
tion and load data under various flight conditions. 
These conditions may best be described with aid of 
Chart 1. 
condition. 
fore omitted). 

A number of variables were recorded on an oscillo- 


The branches represent changes to the basic 
(Condition 5 was unobtainable and is there- 


graph. They were: 
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Calculated rolling pull-out responses; effect of simplifica 


tions 


(1) Aileron position 
(2) Elevator position 
(3) Rudder position 
(4) Normal acceleration 
(5) Fin shear 
(6) Fin bending moment 
(7) Sideslip angle 
(S) Rolling velocity 
(9) Yawing velocity 
(10) Right stabilizer shear 
(11) Right stabilizer bending moment 
(12) Left stabilizer shear 
(13) Left stabilizer bending moment 


In addition, the following items were recorded on the 


photo-observer by a 35-mm. movie camera 


(1) Airspeed (ship’s system 

(2) Altitude 

(3) Outside air temperature 
(4) Fuel quantity remaining 
(5) Airspeed (Prandtl tube 

Run number 


(7) Frame number 
(S) Bank angle 


(9) Time 


It is not proposed to describe the instrumentation de 
tails here. For the interested reader, references 2 and 15 
describe recording instrumentation design procedures 
and the methods used to measure tail loads. 

The flight test program resulted in satisfactory data 
for 5S pertinent rolling pull-out maneuvers which are 


presented in reference 2. These maneuvers are: six of 
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the basic maneuver, four of each of the alternatives, and 
eight with unusual characteristics. The latter were 
obtained during the “build-up” portion of the program, 
a procedure necessary since solutions to the complete 
equations of motion were not available at the time the 
flight program was scheduled to begin. In order to as- 
sure pilot-airplane safety, it was necessary to predict 
the loads to be sustained on an approximate basis, then 
fly several ‘build-up’? maneuvers (increasing g and roll 
rate) until the values specified in the flight program 
were attained. 

At the conclusion of the flight test program, sufficient 
data were at hand to enter what was considered to be 
the most important phase, that of comparing flight 


test with theory. 


COMPARISON OF THEORY AND EXPERIMENT 


Having derived a set of motion equations thought to 
describe the rolling pull-out maneuver adequately, and 
having performed a number of actual rolling pull-out 
flight tests, the next step in the program was to attempt 
to validate the equations of motion through comparison. 
In addition, the shortcomings of equations, at least in 
terms of those quantities which could be readily varied, 
would become evident and could be investigated. 

The set of equations used was essentially that of the 
first simplification of the theoretical part of the pro- 
gram. Due to the fact that during flight, rudder mo- 
tion was inadvertently applied (and recorded), such 
rudder input terms had to be added to the basic rudder- 
locked equations. These were Cy, 5, in the SY equa- 


tion, Ci, 6; in the YZ equation, and Cn, 9r in the >N 
7 Oy 
equation. To further the exactness of the comparison, 


aerodynamic and mass parameters, substituted in the 
equations of motion, were, where appropriate, based on 
the exact conditions of flight test. This was done for 
the initial straight and level lift coefficient, C,, absolute 
angle of attack, a,,, altitude, moments of inertia, and 


tail length. The derivative C,, was calculated for the 


flight c.g. position. For the case where tip tanks were 
installed during the maneuvers, the expressions which 
» Sigs Came Cas G 


describe the derivatives C,,, Ci,, Cn, 


and C., were recalculated to include the “body” con- 


tribution of the tanks. 

Inputs to the equations of motion were those recorded 
for each flight test maneuver selected for comparison 
The conditions selected included all those previously 
noted as test flown, with the exception of conditions 7 
and 11, which appeared to be the least critical. 

On the basis of preliminary comparisons, a further ad- 
justment was made to the derivatives Ci, and C ; to 


r ; 
account for aeroelastic effects, data for which was un- 
available at the time of the initial theoretical studies. 
This essentially amounted to matching the roll rates, 
and adjusting C;, , and matching the initial yaw rates, 


and adjusting C, . All further comparisons and param- 


eter variations were made from this reference. 

As in the case of the initial theoretical calculations, 
preliminary comparison calculations were made using 
digital facilities. The majority of the calculations were 
made using the analog facilities at WADC, which were 
quite adequate for the simplified five degree of freedom 
equations. 

The results of these calculations permitted compar- 
sons, under the flight condition variations previously 
noted, of calculated and measured normal acceleration, 
sideslip angle, vertical and horizontal tail loads, roll 
rate, and yawrate. Although the intent was not to vary 
any of the parameters, a limited additional series of cal- 
culations was carried out in the interest of improving 
agreement between flight test and theory. The param- 
eters varied included C,,,, Cs» Cras Ons I,,, and I,,,. 

This paper presents some of the more interesting re- 
sults of these investigations. The entire comparison 
phase has been reported in reference 3. 
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Fic. 4. Comparison of flight test and calculation; condition 1 
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Fic. 5. Relative effect of parameters varied on maximum side 
slip angle; condition 1. 


Fig. 4 shows the comparison between calculation and 
flight condition 1. Except for the differences already 
noted, this is the basic condition used in the initial 
theoretical study and the basic condition of the flight 
program. It may be noted that the calculated incre- 
mental normal acceleration is somewhat higher than 
flight test. 
somewhat higher. 


Both 8 and the vertical tail load are also 
If normal acceleration agreement is 
obtained (in this case it was done by changing Cm)» 
then 8 and .Vy amplitudes show reasonable agreement. 
The second peak in the Ny response is much larger, 
suggesting that the rudder effectiveness might be too 
large. No attempt was made to validate this. The 
extent of parameter variation attempted is illustrated 
in Fig. 5, where the effects are presented in terms of the 
change in maximum sideslip angle, as a function of the 
change in parameter value. A strong coupling effect is 
evidenced by the large variation of maximum sideslip 
angle with Cn, The obvious effect of Cr, is noticed 


as is the extreme importance of Ci, . Thus, two of the 
°a 


derivatives, (C,, and C,,. ), which exhibit nonlinearities 
Qa oa 
and which were comparatively difficult to estimate, 
show a very strong effect on the motion. 
Fig. 6 shows the comparison for condition 12; similar 
to condition 1 except that instead of the elevator lead- 


MANEUVERS 5138 


ing the aileron, the two control inputs occur somewhat 
simultaneously. Generally good agreement was ob- 
tained for this condition, although adjustment of An, 
for agreement at the first peak would result in a slightly 
greater unconservative sideslip response than shown 
here. No parameter changes were investigated for this 
case. 

Condition 13, again similar to condition | except for 
the control inputs, which in this case includes rudder 
application in an attempt to regain coordination, is 
compared to calculation in Fig. 7. Agreement again is 
reasonable except in the case of yawing velocity where, 
perhaps because of too large a C,,., the calculated r is 


much larger than flight. Reduction of C,. would un 


doubtedly reduce the maximum sideslip angle and ver 
tical tail loads, thus producing slightly unconservative 
values for these variables. 

In all cases the flight test forward speed change 
is shown alone, since this was neglected in these calcu 
lations. 

As was previously implied, no attempt was made in 
this study to investigate the effects of all important 
derivatives and mass parameters separately or in all 
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Fic. 6. Comparison of flight test and calculation; condition 12 
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reasonable combinations which may have produced 
near perfect agreement with flight test. Such effort 
as was directed along these lines was done so on the 
basis of those variables which were thought to be least 
accurately estimated. 

Important results were obtained from this compart- 
son and limited parameter variation study. 
crepancies noted between flight test and calculation 
were similar to those encountered in small disturbance 
maneuver comparisons (basically attributed to inade 
quate aerodynamic data), thus lending credence to the 
five degree of freedom method of calculation proposed, 
but indicating that, because of the importance of aero 
dynamic coupling, it might be necessary to estimate 
aerodynamic data with greater accuracy and precision. 
The initial theoretical study pointed to the need of 
five degrees of freedom to describe large disturbance 
motion. This was further emphasized in the param 
eter study by the strong effect of longitudinal param 
eters on lateral motion. 

The promising results of the overall initial investiga- 
tion were to provide a firm background for an extended 
program designed to investigate the “less conventional’ 
airframe, where the effects of more complex aerody- 
namic coupling and, in particular, the effects of more 
extreme inertia coupling could be studied. 


A SIMPLE APPROXIMATION 


While work was proceeding on the comparison phase 
of the program, some thought was being given to the 
derivation of a relatively simple expression for maximum 
sideslip angle which would include the essential param- 
eters, yet would dispense with the arduous computa- 
tion associated with the five degree of freedom set of 
equations. 

Basically, it was felt that several coupling terms 
could be included without undue complexity. Most 
— [,)pq, appearing 
This term includes 


important was the coupling term (/, 
in the yawing moment equation. 
three important rolling pull-out parameters, one of 
which, of particular interest, is the moment of inertia 
difference. Emphasis was placed on this difference, 
since the trend in modern fighter design indicated larger 
and larger values for it. 

The approximation should be defined in terms of basic 
parameters which could be readily assumed for pre- 
liminary prediction of maximum sideslip angle. Load 
factor, roll rate, and basic flight condition appeared to 
be most promising. 

The yawing moment equation was reduced to the 


following: 


ce 


Y 


— I,)pq = gSbC,, 8 + gSbC,,,(a)(pb/2U) + 
gSbC,,. (a) 6, 
oq 
Letting the yawing moments due to p and 6, (and a) 
be represented by AC,(a), 
(I, — Iz)pq/qSb = C,8 + AC,(a) 


g is approximated in the following manner. The nor- 
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mal force equation in terms of load factor may | 


proximated as 
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—g + Bp = (g/U)An g/l) ) 
or a 
q = 3p _ g [ i 1) 
Substituting in the yawing moment equation, 
Ud, — I,)p[Bp + (g/U)(n -— 1 
- C,.8 + AC, (a 
good ” 
Ud, — 1,)8p? + Ud, — 1,)(g/U)(n b I 
qb 
C4,0 + AC de 
(I, — I,)8p 
en - U, 6 = 1 
gSb J rh 
; f. [.)(g/U)(n | wes 
At Ae = of a 
qs0 
Stal 
AC, (a) — [U, — 1,)(¢/U)(n 1) p/qSb 
(1, — I,)/qSb p* — c. 
= > 
: ; In 
If steady-state values for » and p are assumed, along 
with a dynamic overshoot factor of 2 (assumption oj oe 
. . . . tT 
zero damping), the maximum sideslip angle becomes 7 
con: 
" Ud, — I,)(g (mss — 1)p It 
AC,(a) — : 
> gSb oes 
Bmar = 2 . 
(UZ, —_, gb |p, 7 — ( , 
A word must be said about the AC,,(a@) term. As pre Thi: 
viously noted, 
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AC, (a) = C,,,(a@) (pb 2U) + Cy, (aba 


and C,,. angle of attack functions must 


The actual C ; 
be determined for the individual airplane, and the angle 
of attack chosen must be commensurate with the steady 
state load factor chosen—4.e., 


\(W/qS) 


“a 

In addition, in the case of C,,. (@)6,, a 6, must be chosen 

commensurate with the roll rate, p. Such may be de- 

termined by calculation for the flight condition under 
consideration. 

It may be noted from the expression for 8,,,, that 8 


goes to infinity when 


I, — 1,)/gSb]Pu? = Cn, 
This may be caused by a combination of large (/, — /,)» 


, and low C,,,—all characteristics of modern fighter 
For a given airplane there exists a p,, for 


This has been termed the 


airplanes. 
which the above occurs. 
critical roll rate, and more will be noted about this con- 
dition later, on the basis of more recent work than that 
described above. 

Very little work was done with this expression at the 
time of its derivation, except to note the effect of mo- 
ment of inertia difference in comparison with those ap- 
proximations which did not include this effect. It is 
ipparent that the approximation is valid at least within 
respectable limits for what was then considered a con- 
ventional airplane. Fig. 8 shows a comparison of the 
approximation with the exact solution as described in a 
previous section, and with one of the flight conditions. 
Application to a swept-wing modern fighter with char- 
acteristics noted above indicates limitations which for 
approximation. 


practical invalidates the 


However, in the case of other types of modern aircraft, 


purposes 


such as a light or medium tactical bomber, possessing 
characteristics similar to the earlier fighter types, and 
required because of mission to perform large-disturbance 
maneuvers, the application of this approximation may 


be justifiable. 


STuDY OF A HIGH-SPEED AIRPLANE 


"he modern fighter is a different airplane from the F- 
SUA for which the five degree of freedom equations were 


DISTURBANCE MANEUVERS O15 


developed. Present-day trends are toward low aspect 
ratio wings, high wing loadings, and mass concentrated 
in the fuselage. Accordingly, the dynamics of a typical 
swept-wing fighter were studied in general rolling ma- 
neuvers, emphasizing variations in wing plan form, ver 
tical tail size, mass distribution, flight condition, aero 
dynamic derivatives, and control inputs. The re 
mainder of this paper is concerned with the results of 
this study.‘ 

The configuration chosen for study is compared in 
Fig. 9 with the F-SOA and was based on a prototype 
version of a typical swept-wing fighter. Contemporary 
trends in geometry and mass characteristics are readily 
apparent. The aerodynamic derivatives were esti- 
mated for this configuration and include power effects 
and aeroelasticity. The results are presented in Fig. 
10 for JJ = 0.9 and 20,000 ft. altitude. These deriva 
tives are referred to stability axes. Only large non 
linearities are included, but no attempt is made to 
study the region of the stall. The variations of C,., 
Ci», Cyp, and C;,, with @ are typical of airplanes with 
low aspect ratio swept wings. The forces and moments 
due to aileron deflection are peculiar to the particular 
lateral control system used. The equations of motion 


are written with respect to body axes as 


mV(B — ap + r) = qSC, + mg cos@sind (6a 
mVla + Bp — gq) = gSC + mg cos cos @ (6b 


= gS0C 6c 


L.pt+ Ud, — I,er — 12:(¢ + pq 


(p>? —r : gs¢ c. 6d 


Lg+ (UI, —I)pr+! 


Lt + (I, — 1) pq — Lp — gr) = qSbC, (Ge 


These are Eqs. (1) with the following simplifications: 


= w/u, B=v/V, amv constant 


where u is the total component of velocity along the x 
axis. Airplane orientation in space (yy, @, and @) is 
given by Eqs. (3). The incremental reading from level 
flight of a normal accelerometer placed at the center 
of gravity of the airplane is given by 

—An, = —(V/g)(a& + Bp — g) + cos@cos@ — 1 (7 


where the negative, —An,, is used so that a positive 
‘“‘g’’ will correspond te a pull-up. In level flight 
—An, = 0. 


The equations may be conveniently written as 
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Fic. 11. Comparison of derivatives in stability and body-axis 
systems. 

B-—ap+r=ye (g/V) cos@sin o (Sa) 

a+ Bb —gq =2+ (g/V) cos 86 cos @ (Sb) 


pt (Uz — [,)/Trlar — U22/Tr)(* + pq) = 1 (8e) 
g+ (U, — 1,)/I,]pr + U22/T,)(p? — 7?) = m= (8d) 


*+ (UZ, — 1[,)/l.|bq — U2:/12)(p — gr) =n (8e) 


where 
y Y/mV = (pSV/2m)C, 
gs = Z/mV = (pSV/2m)C, 
Ll = L/I, = (qS0/I,)C, 


M/TI, = (qSc hCG 
n = N/I, = (qSb/I,)C, 


It is necessary to transform the aerodynamic forces 
and moments from stability axes, the system in which 
they are usually given, to body axes, the system in 
which the equations of motion are written. This sub- 
ject, the transformation of data from one axis system 
to another, is a very important one and is discussed 
extensively in reference 4. The equations for trans- 
forming the derivatives are outlined in Appendix 1. 
In general, the derivatives in the stability axis system 
may be radically different from the corresponding 
derivatives in the body axis system. The derivatives 
which show the greatest difference are usually the yaw- 
ing moment derivatives, and particularly Cry» Cue 


and C,, which are important ones. In order to illus- 
°aq 


trate this point the latter two derivatives are shown 
in Fig. 11 in both axis systems for the swept-wing air- 
plane at 1/ = 0.9 and 20,000 ft. altitude. 

For the study airplane, the body axes are chosen so 
that they coincide with the stability axes for straight 
and level flight at the reference flight condition. After 
transforming the aerodynamic data to body axes, the 
aerodynamic forces and moments, including engine 
gyroscopic moments, may be expressed in nondimen- 


sional form as 


y= y,(a)p + yer + V8 + Vsa5e (9a) 


Z = 2G + 2G + Bq2(a)bg” + 25-0 (9b) 


l=1,p + liq + 1(a)r + 13(a@)B + Isa(a, 6a)6. (9C) 
m= mM,p + mg + mr + maga + Ma + 


Miq2?(@)6q" + ™5,6¢ (9d) 
n= Nn,(a)p + mq + nr + ngB + Ngq(a, 5a)b. (Ye) 
where the variable coefficients for the study airplane 
resulting from nonlinear aerodynamic derivatives are 
expressed in polynomials as 
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Eqs. (8), (9), and (10) comprise the final five degree oi 
freedom equations of motion for large disturbances 
from level flight, and the required kinematic relations 
for 6 and @ are given by Eqs. (3). These are the 
equations used to compute F-SOA solutions for com 
parison with flight test, except that the small angle ap- 
proximations, sin 6 = @ and cos @ = 1, were included 
and the nonlinear derivatives used for the F-SOA were 
l(a), mM,(a@), Np(a), n-(a), and n;,(@). 

Three flight conditions were studies for the swept 
wing fighter: 


Condition A: M = 0.9, h = 20,000 ft. 
Condition B: M = 0.7, sea level 
Condition C: M = 0.7, # = 30,000 ft. 


and the coefficients in the equations of motion are give! 
in Appendix 2. 

Analog solutions for condition A are presented in Fig 
12, comparing an aileron roll with a rolling pull-out 
initiated from a 5.11g pull-up (maximum unsymmetrical 
Aileron deflections are selected to givé 
The aileron roll solu- 


load factor). 
comparable maximum roll rates. 
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hibits extreme sideslip angles and noticeable coupling 
in the a and — An, records. 

Solutions for condition B are presented in Fig. 13 
comparing two rolling pull-outs initiated from a 5.11g 
pull-up. These solutions are even more extreme than 
those for condition A and show little resemblance to 
what one normally expects in a rolling maneuver, par- 
ticularly for the 30° aileron deflection. The sideslip, 
8, oscillates about a steadily growing value and even- 
tually reaches a very large steady state. Large ampli- 
tude longitudinal motions are excited, normal accelera- 
tion — An, varying between 0 and 8 g’s, a severe ma- 
neuver by any standard. Comparing the two solutions 
in Fig. 13, we see that nonlinearity is the predominant 
characteristic. 

Solutions for condition C are presented in Fig. 14 
comparing two rolling pull-outs initiated from a 2g 
pull-up. The maximum load factor is limited by the 
angle of attack in this flight condition. The solution 
with 20° aileron deflection appears similar to the roll- 
ing pull-out solution shown in Fig. 12 for condition A, 
but the solution with 30° aileron deflection seems com- 
pletely divorced from a rational maneuver. Just 2.5 
sec. after the ailerons are deflected, the roll rate is 
—340° /sec., the angle of attack is —15°, normal ac- 
celeration (— An.) is —3.5 g’s, and the sideslip is 5 


having previously exceeded —10 Comparing the 
two solutions, we find that for 6, = 20° the average 
roll rate is about —60°/sec., whereas for 6, = 30° the 


average roll rate is about — 300° /sec. 

The solutions presented in Figs. 12, 13, and 14 are 
typical of those obtained at the beginning of this inves- 
tigation. The results were so different from anything 
previously encountered that their validity was ques- 


perienced severe maneuvers at extreme altitudes, but 
the study airplane was a reasonable model with the 
external geometry of a current prototype fighter. 
However, investigation showed that the prototype had 
experienced severe uncontrollable motions in flight 
and that several other prototype airplanes were experi 
encing similar difficulties (see reference 9 A flight 
test record of a rolling maneuver was obtained for the 
prototype fighter. The appropriate inertial character 
istics and flight condition were used but with the aero 
dynamic derivatives of this study to obtain the con 
stants in the equations of motion. Satisfactory agree- 
ment between the flight record and the analog solution 
was obtained using a rough approximation for the con 
trol inputs recorded in flight. This shows that the 
unexpected nature of these computed responses is in 
fact reasonable. The question becomes, why do they 


occur ° 


STUDY OF LARGE AMPLITUDE RESPONSI 


One method of obtaining an understanding of the 
solution of a set of differential equations is to examine 
the contributions of the individual terms. The rolling 
pull-out solution of Fig. 12 is substituted in Eqs. (5S), 
rearranged with only the accelerations 8, a, p, g, and 7 
on the left-hand side. The various terms in each 
equation, the inertial terms, the gravitational terms, and 
the lumped aerodynamic terms, are plotted in Fig. 15. 
The product of inertia terms are negligible in this solu 
tion because the body axes nearly coincide with the 
principal axes, so that n > a and /,, > 0. The 57.5 
factor is necessary because the results are presented in 
units of degrees. In order to complete the picture the 
components of the aerodynamic forces and moments, 


as given in Eqs. (9) and (10), are shown in Fig. 16 
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lerms not presented are of an order of magnitude 


smaller than the ones that are shown. 

From the 8 equation we see that the force terms are 
small and that 6 is dominated by the ap and —r 
terms. Similarly, @ is dominated by the —p and q 
terms. Since these four terms are essentially kine- 
matic and just account for the rotation of the body 
axes, then the motion seems primarily one of rotation 
defined by the three rotational degrees of freedom. The 
roll rate is primarily a function of the aerodynamic roll- 
ing moments as seen from the p equation. However, 
the pitching and yawing velocities are heavily influenced 
by the —[(/, — /,)/I,|pr and —[(U, — I,)/I:]|pq terms 
in the g and / equations, respectively. These two terms 
have been called the gyroscopic terms and are the real 
villains. ware at the components of 7, we see that 
for 2.0 <t< 3.0 sec. the —[(/, — /;) 
exactly cancels the yawing moments, 1. 
16, we see that 7 is largely 


I,|pg term almost 
Examining 
the components of n, in Fig. 
n,B8 since n,(a)p and n,;,(a, 6,)6, about cancel. Thus, 
the gyroscopic term in the yawing moment equation 
tends to cancel the directional stability. Since the 
“net yawing moments’ are nearly zero, the sideslip 
builds up to large values because of the ap term in the 
B equation. Animportant point is brought out. Since 
the body axes for this solution are principal axes, the a 
in the ap term represents the angle of attack of the 

principal axis. Thus, the location of the principal 
axes with respect to the flight path has a direct bearing 
on the magnitude of the response in a rolling maneuver. 
If the principal axis coincides with the flight path 
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(yn = 0), this mode of generating sideslip does not exist 
This is one reason why the sideslip is comparatively 
small in the aileron-roll solution in Fig. 12. 

The results shown in Figs. 15 and 16 suggest that 
many terms in the equations of motion can be neglected 
without greatly affecting the solution. A simplified 
set of equations of motion is derived on the following 
bases: (1) all aerodynamic terms not shown in Fig. 16 
are neglected, (2) the gravitational forces are neglected, 
and (3) the nonlinearity of rolling and yawing moments 
with aileron deflection is neglected. The equations 
are written with respect to principal axes as 


B-—ap+r = y,0 lla 
a+ pp — ¢ = £0 L1b 


(J. — [,)/Ielar = lpp + 1g(a@)B + lsq(a)bq (Ile 


(1, + 1,)/I,|pr = mq + maa + m,;6, (11d 


[(I, — Ir)/Iz\pq 


n,¥ + NgB + 
N,»(a)p + N5,( a) 6q lle 


The damping in pitch and yaw are retained, even though 
their normal effect is small, since they are often the 
source of stability augmentation. This simplification 
of the equations is treated with considerable detail in 
reference 4, but these equations should give relativel) 
good quantitative answers if the following conditions 
are met: (1) principal axes are used as the reference, 
otherwise product of inertia terms must be included 
(2) the airspeed is high enough so that the wing angle ol 
attack in level flight is small; and (3) /;, and m,, ate 


wl 


th 
th 





noments 


ot exist 


iratively 


est that 
eglected 
mplified 
owing 
Fig. 16 
elected, 
10Oments 
juations 


(Lle 


though 
ten the 
fication 
etail in 
atively 
ditions 
erence, 
luded; 
negle of 


N5q are 


AIRPLANE DYNAMICS IN LARGE 


idjusted for the usual decrease in aileron effectiveness 
at large deflections. Obviously, if the magnitude of 
the motion exceeds the valid range for the aerodynamic 
representation, go¢ | results cannot be expected. Fig. 
17 compares a solution computed using the full equa 
tions with one computed using the simplified equations. 
[his is the rolling pull-out solution shown in Fig. 12 
ind the one used to prepare Figs. 15 and 16. These 
simplified equations will usually be more accurate than 
the numbers used in them. Not only are they com 
paratively easy to solve, requiring a minimum of 
numerical data and analog equipment, but they may 


be profitably studied using analytical techniques. 


\NALYSIS OF THE EFFECT OF ROLL RATE 


Analysis of the equations of motion for steady rolling 
explains much of the behavior of airplanes in rolling 
maneuvers, and several references are available on the 
subject This approach was first used, to the au- 
thors’ knowledge, by White, Uddenberg, Murray, and 
Graham’ who considered the stability, the frequency 
response, and the steady-state solution of the equations. 
A two degree of freedom approach (freedom in pitch 
and vaw) was developed by Phillips,® and Bergrun and 
Nickel’ verified Phillips’ results using a freely falling 
model (body-tail combination) with a roll-stabilization 
system. Phillips considered only the stability (linear- 
ized equations in pitch and yaw) for steady rolling 
about the flight path, but the analysis is quite complete. 
A comprehensive theory for the dynamics of rolling 
missiles was developed by Nicolaides* for bisymmetric 
,=, ¢ = ¢, - &. * 


a 3 a 


configurations —41.e. 


b/c)C,,, ete. Nicolaides uses a complex variable 


3 
notation-—particularly applicable to bisymmetric air- 
craft--which presents the effects of rolling in a particu- 
larly clear manner, and his paper includes both steady 
rolling and roll dynamics. Since the effects of steady 
rolling are adequately covered in references 5-7, we 
will only include the important results as needed. 

The first step is to put the equations of motion in a 
more meaningful form by eliminating some of the vari- 
ables. If we solve Eqs. (lla) and (11b) for 7 and gq, 
differentiate to obtain * and g, then r and g can be 
eliminated from the pitch and yaw equations, (11d) 


and (lle), to yield 


&+ (—2, — m,)a + (—m, + 2M, + Op*)a + 
| — Q)ppB + (Qyg — m,)p8 + p8 = m,,6. (12a) 
B + Yg — n,)B + (ng + ygn, — Rp*)B — 


(1+ R)pa + (R.. + n,)pa — pat+ny, 


a)p = 
—n;,(a)6, (12b) 

where 

QO = (J, — I,)/I, R= (,.— £)/I, 

lf we constrain the motion to one of steady rolling about 

the x principal axis, with p; = steady roll rate, then 

the equations can be linearized for small disturbances 

in the remaining four degrees of freedom. If we as- 
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fied and full equations of motion 


sume that ,(a@)p and n;,(a@)é, are the inputs to the 
yawing moment equation, then the pitch and yaw 


equations may be written as 


a+ b,a4+ k,a + (1 — Q)pis 4 


(Ovg — m,)pi38 =m, (Lda) 


8 i 538 . k 38 = (1 _ R Pia T 
(Re, + n,)pia -n, (13b) 


where the following notation has been introduced 


b, = -2,—m 
k, = m, + 2m, + Qp 
03 = —V¥g— Nn 
kz; = ng t+ yen, — Rp 
mM. = m;,6 
No = Nn,(a)p + nz,(a@)6 
We notice that, for p; = 0, b, and k, are the usual lon- 


gitudinal short-period total damping and _ stiffness 
(usually includes m, which is neglected in these equa- 
tions), and 6; and kg are the usual approximations for 
the Dutch-roll damping and stiffness. 

Phillips, in reference 6, describes the stability of the 
system represented by Eqs. (13a) and (13b) in terms 
of the characteristic equation. Principal results are 
as follows. (1) For p; = O there is one frequency of 
oscillation in pitch and one frequency of oscillation in 
yaw. (2) For steady rolling, there are two frequencies 


of oscillation in pitch and the same two frequencies of 
oscillation in yaw. (3) One of these frequencies is 
higher than either of the nonrolling frequencies, and the 
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Fic. 18. Stability parameters in equations for steady rolling 


other is lower than either of the nonrolling frequencies. 
(4) As the roll rate approaches zero, the high fre- 
quency approaches the higher nonrolling frequency, and 
the low frequency approaches the lower nonrolling 
frequency. 

If the linear differential Eqs. (13a) and (13b) are 
solved for their characteristic equation, it is found to 
be of fourth degree. The constant term in this quartic, 
which equals the product of the undamped natural fre- 
quencies squared—1.e., the product of the roots—is 
designated by / and may be written as 


E = kukg — (Qyg — m,)(R2, + mr) p1° (14) 


a 


and the system is unstable if HE < 0. The values of 
k,, kg, and E are plotted versus roll rate in Fig. 18 for 
the swept-wing airplane in condition A. For 3 < pi < 
4.8 rad./sec. (1.e., 170 < pi < 275 deg./sec.) the system 
is unstable. The roll rates for which & = 0 are called 
the critical roll rates (p,,;,). 

The following simple criteria is proposed by Phillips. 


Since 
E = kikg 
k, = —mM,, + Op;* 
kg = ng — Rp 


then for stability 


—m, > —Qpr (15a) 


and n3 > Rp, 15b 


If these are reduced to more common notation, then th 
following criteria are obtained: 


b< VW (1/2)pV7ScC,,,/U. — 12) = Derts Iba 
p< V (1/2)pV?S0C,./ Ty, — Lr) & Perit l6b 


Usually Eq. (16b) defines the lower critical roll rat 
These can be considered as minimum criteria for C 
and Cus or they can be considered as criteria for th 
allowable roll rate. Zimmerman, in reference 9, dis. 
cusses the trends in fighter development that have led 
to this problem. The trends may be summarized as 
higher wing-loadings, higher altitudes, and weight con- 
centrated more in the fuselage. 

The analysis of stability in pitch and yaw for steady 
rolling results in attractively simple expressions and 
criteria, but gives no knowledge of the amplitudes oj 
motion in practical rolling maneuvers. More meaning- 
ful results can be obtained by considering the steady- 
state solution of the equations of motion—.e., &@ = 
B=g=r=p=0. 

The steady-state solutions for a and 6 can be obtained 


from Eq. (13) as 
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Fic. 19. Steady-state solution of equations of motion. 
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AIRPLANE DYNAMICS IN LARGE 


kgm. + (Qy3g — m,)pin _ 
y = ; (l7a 
kks — (Qvg — m,)(R2, + n,)pr 
—kn,. — (Rz, + n-)pim 7 
= (17b) 


kks — (Ovg — m,)(Rz, + n,)pr 


The vawing and pitching velocities are obtained from 


Eqs. (lla and (11lb) as 
ry = AssP~i + V30 (1l7e¢ 
q = p—, Pi = 2a (17d) 


The rolling equation may be written, with the gr term 


neglected since it is usually small, as 
0 = l pp + l3(a)3 + | (17e) 


We note that the denominator in Eqs. (17a) and (17b) 
is the constant term in the characteristic equation, £, 
of the steady-rolling analysis. 

A complete solution of these steady-state equations, 
including the nonlinear derivatives ”,(a@) and n,,(a 
in m, and /3;(@), requires a lengthy iterative procedure. 
For illustrative purposes the steady state is computed 


using 

mM. = Msdbe, a, = 0, 1. = 15.60 
This allows the steady state to be presented as a,,/ —4, 
and 8,,/—6,, though 6, is still not proportional to 6, 


for any given value of p;. The results are presented 
in Fig. 19 for the swept-wing airplane in condition A. 

The instability predicted by the linearized equations 
for steady rolling (£ < 0 as shown in Fig. 18) appears 
in the full equations as a new equilibrium state, charac- 
terized by large a,, and @,,. Discontinuities are found 
at two roll rates corresponding to & = OQ, and are the 
critical roll rates. The steady-state solution is inde- 
terminate at the critical roll rates and the dynamic 
solution tends to “jump” from one branch to the next. 

A number of important aspects of the dynamic re- 
sponse are explained by the steady-state solution. In 
general, the origin of the large amplitude motions en- 
countered at high roll rates is seen to be primarily the 
nonlinear inertial terms in the pitch and yaw equations. 
We see from Fig. 18 that if k, = 0 and kz = 0 at the 
same roll rate, then the amplitude of the steady state is 
minimized, though it still may be large. This corre- 
sponds to the condition that 


(—cCng)/(Iz — I) & bC,,/(I, — 12) 


or, approximating further, that the natural frequencies 
in pitch and yaw are equal. In addition, Fig. 18 shows 
that increasing the overall damping through any of the 
derivatives m,, ,, 54, OF Y3 Will decrease the steady- 
state response, since the increment between E and k,k, 
is thus increased [see Eq. (17) ]. 

The aspects dealt with so far concern factors in the 
denominator of Eqs. (17a) and (17b). The numerators 
of these two equations are comprised of two types of 
terms: (1) the pitching moment due to longitudinal 


control, and (2) the yawing moments due to roll and 
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due to ailerons (also the rudder control, which is neg 
lected in this study). The longitudinal control in ef 
fect defines the angle of attack of the x principal axis 
for the nonroliing airplane, and the adverse yaw deriva 
tives are known to be very sensitive to the location of 
the reference axes (principal axes in these equations). 
Thus the angle of attack of the principal axes and the 
adverse yaw derivatives with respect to these axes are 
determining factors in the response. 

The steady-state relationship between roll rate and 
aileron deflection as shown in Fig. 19 demonstrates one 
of the most significant aspects of the problem 
namely, that the roll rate is not uniquely defined by the 
aileron deflection. In fact, four possible autorotational 
roll rates exist—1i.e., conditions of steady rolling with 
no aileron deflection. From Eq. (l7e) it is apparent 
that /;(@) is responsible for this behavior. No exten 
sive study has been made of these autorotational con 
ditions, but it has been noticed that once the lower 
critical roll rate is exceeded, there is a marked tendency 
for the roll rate to jump to the higher critical roll rate. 
This indicates that the lower autorotational roll rates 
are unstable whereas at least one of the higher roll rates 
may be stable. The important point is that aileron 
control power does not limit the roll rate and that a 
dynamic solution including both inertial and aerody- 
namic nonlinearities is required to estimate the roll 
rates in a rolling maneuver. 

Stability in steady-rolling flight and the steady-state 
solution of the simplified equations {Eq. (11) ] are seen 
to be concepts useful in understanding the nature of 
the equations of motion and their solution. Though 
their application to quantitative calculation of dynamic 
responses is limited, these two concepts are good tools 
for studying the qualitative effects of various param- 
eters. However, the steady-rolling equations and 
the steady-state solution can yield quantitative answers 
for some problems. 

Referring back to Fig. S, we see that Eq. (5) pre- 
dicted 8,,,, in a rolling maneuver for the F-SOA air- 
plane with creditable accuracy. 8,.,, as given by Eq. 
(5), is just twice 8,, as given by Eq. (17b) with m, = 
n, = V0 = 0. 
expected to give good results if (1) the two critical roll 
i.e., if the natural frequency 


The simple expression for 8, can be 


rates are well separated 
in pitch is much higher than the natural frequency in 
yaw—and (2) if the roll rates are well below critical. 
These conditions will be met by most airplanes having 
straight wings with high aspect ratios. 

The equations of motion written for small disturb- 
ances from steady rolling may also be profitably used 
in the design of automatic stabilization and flight con- 
trol systems. Since presumably the stabilization sys- 
tem will restrict the motion to small amplitudes, then 
the linearized equations should be valid. The root- 
locus technique can be used to plot the loci of the 
zeroes and poles with roll rate as a parameter. The 
stabilization system must then provide satisfactory 
response throughout the range of applicable roll rates 
Roll dynamics may be included by allowing small per 
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Effect of varying the location of the principal axes on 
roll response 


Fic. 21 


turbations in the roll rate. This linear approach does 
not provide means for computing actual responses in a 
given maneuver. Only the full equations of motion 


can provide such responses. 


ANALOG SOLUTIONS 


Means for understanding and interpreting the dy- 
namic response in rolling maneuvers have been de- 
veloped in the previous section. In this section we will 
examine the effects of various parameters on the com- 
puted responses for the swept-wing airplane in flight 
condition A. The full equations of motion |Eqs. (S) | 
have been used to compute these responses. 

The source of the large-amplitude motions in rolling 
maneuvers has been traced to the nonlinear inertial 
terms (inertial coupling) in the pitch and yaw equa- 


tions: 


(J, — I.)/T,|pr 
[Z, — 1,)/T2]pq 


in the pitch equation 
in the yaw equation 
The effect of increasing /, by a factor of three is shown 


in Fig. 20. Coupling of pitch and yaw is greatly re- 


ICAL SCIENCES FUL 
duced and the response is similar to what would be ex 
pected from linear equations. 

The location of the principal axes has been stressed 
as important. The icates 
that n (angle of attack of the x principal axis) acts like 


steady-state solution ind 


f 44 


an input in a roll and dominates the magnitude of th 


response. <Aileron rolls for three locations of th 
principal axes are compared in Fig. 21 Phese varia 
tions, +3° and 5° from the bas urplane, are not 


large, and their effect on the response is nothing short 
of startling. In all three cases the roll rate is near the 
critically roll rates, but for the basic airplane there is 
no mechanism for exciting large amplitude motions. 


The rolling responses for a series of initial pull-ups ar 


summarized in Fig. 22 by giving the extreme values 
(after the first maximum) of a, 3, and p obtained in 
each solution. The reference values of @ and p 

1.€., Qrer and P,, are essentially the angle of attack 


and roll rate that would be obtained from the uncoupled 
lateral and longitudinal equations. The largest excur- 
sions in a and @ are seen to correspond with roll rates 
near critical, large @ occurring for p,., near the low 


Perit (4 g's), and large @ occurring for p,,, near the high 


Perir (O g's) and where the roll rate oscillates about the 
high critical roll rate (3 g's). These results correlate 
with the steady-state results shown in Fig. 19. The ex- 
treme roll rates obtained for 3g reflect the autorotational 
tendency. 

Variations in wing and tail geometry have been stud- 
ied and the most significant parameter was found to be 


the vertical tail size. Results for this variation ar 
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The influence of the inertia-coupling terms may be 
shown particularly well if the yawing moment due to oes ee Ff Ts 
= ae a re ‘ t ” e t _ = 
roll is favorable and large. In Figs. 26 three solutions 
are compared for the swept-wing airplane but with a 
positive C,, as in Fig. 25. In the first solution, the x 
} inertia-coupling terms 
[(U, — I.)/T,|pr [J, — 1,)/I2\bq ; ; 
are set to zero and the response is what one might expect ; a AS ans 
irom a linear airplane. In the second solution the in- ——— —— 
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equations. 


normal acceleration. Control motions are also critical 
In Fig. 27 three relative time lags between 


The eleva- 


quantities. 
the elevator and aileron inputs are shown. 
tor reversal at high roll rates—i.e., change in sign of 
that was noted in the steady-state solution is 


Elevator 


a/6, 
clearly demonstrated by these responses. 
deflections applied during the roll are generally found to 
be critical. This also applies to aileron deflections 
since in some maneuvers the largest sideslip angles 
have been encountered after the ailerons have been 
neutralized to stop the roll. Though control motions 
are important, the airplane's orientation in space is of 
secondary significance. In reference + two solutions 
are compared, one for a roll initiated from a pull-up 
and one initiated from a steady turn, and the two solu- 
tions are identical except for spacial orientation—4.e., 
y, 6, and #. Of course, this result might be expected 
since the simplified (11)], which 
neglect gravity, agree well with the complete equations. 


|Eqs. 


equations 


Basic NATURE OF MOTION 


The basic nature of motions in a rolling maneuver is 
to some extent described by the steady-rolling and 
steady-state analyses developed earlier. However, a 
clear understanding of the nonlinear phenomena and 
the relationship between the aerodynamic and inertial 
forces may be gained by examining the force-free mo- 
tions of a rigid body—e.g., motions in a vacuum. 

For symmetrical bodies (/, = /,), the force-free equa- 
tions can be solved using elementary functions. 

Iq = UI, — I[,)pr 
Ir = —(U, — [r)pq 
[,p = 0 
The variables may be separated by differentiation and 
substitution to yield 
ies § to 
G = — [Uy — I;) I, |ps q 


r = —([(Z, — I)/T,]p$*r 


SCIENCES JULY, 195 


These equations describe simple harmonic motion 


A cos wl, P= sin wl 


q — 
where 


w = (ly —1,)/L)p 
p = 


A = constant 


constant 


The magnitude of the total angular velocity is a cor 
stant but the position of the angular velocity vector 
rotates (precesses) in the body with frequency w. Thus 
steady rotation can occur only about a principal axis 
(4 = 0), and rolling about a nonprincipal axis will in- 
troduce precession, causing excursions in a and 3 
An important point is made—the force-free motion js 
neutrally stable. The body is no more prone to rotate 
about a principal axis than it is to rotate with preces 
sion. 

For unsymmetrical bodies (/, # /.), the equations may 
be solved in terms of elliptic functions, but a geometrical 
construction in three dimensions, known as Poinsot’s 
construction, can be used to show the character of the 
this 
may refer to reference 10 where the equations are dis- 


motion. For the details of Poinsot motion one 
cussed and a pictorial description of Poinsot’s model is 
shown. The principal result is the same as for the 
symmetrical body, steady rotation can exist only about 
a principal axis. For any other condition the body 
will rotate in a predictable but nonuniform fashion, 
causing excursions in @ and 8. If /, is nearly equal to 
/,, aS in many recent airplanes, then the motion will be 
nearly periodic and very similar to the symmetrical 
case. 

The essence of the motion in rolling maneuvers can be 
described by considering two extremes (force-free and 
inertialess). 

(1) For a rolling airplane with aerodynamic forces 


but zero moments of inertia: 
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AIRPLANE DYNAMICS IN 
a = ( Cr) (18a) 
8 = [C,, 6. + (6/2V)C,,p)/Cu, (1S8b) 
») For a rolling airplane with inertia but with zero 
porces 
a = —pp Tq (19a 
6=ap-r 19b 


and a and 8 can be constant only if p is about a princi 
pal axis and g and r are zero. 

~ Thus, for a rolling airplane, if the principal axis coin- 
cides with the flight path, as described by Eq. (1S 
steady rolling is possible without otherwise disturbing 
the airplane. However, for any other condition Eqs. 
IS) and (19) represent a conflict, and the airplane will 
experience lateral and longitudinal disturbances. The 
magnitude of these disturbances will depend primarily 
on the roll rate and the static stabilities. 

Two very simple models may be used to demonstrate 
the inertia-coupling phenomena. The first is a rod, 
pin-connected to a shaft or axis of rotation and re 
strained to the shaft by aspring. The rod represents the 
airplane, the spring represents static stability, and 
the shaft represents the component of rotation along the 
flight path. By rotating the shaft at various speeds 
the existence of a critical roll rate and large-amplitude 
steady states can be shown. The equation of motion 
for this model is very simple and demonstrates stabil- 
ity and steady state as a function of roll rate. 

The second model may be made by hanging a small 
toy airplane (preferably metal and with short wings) 
by its tail with a rubber band. The rubber band and 
gravity supply static stability in pitch and yaw. If 
the model is wound up on the rubber band and released, 
it will demonstrate the inertia-coupling phenomena by 
revealing the different modes of motion associated with 
different roll rates. 

Time and space do not allow a complete presentation 
here of this type of approach which utilizes elemental 
equations. However, recent studies at Cornell Labora- 
tory suggest that further effort along these lines may 


prove very valuable. 


CONCLUSION 


In this paper we have presented a five degree of free- 
dom approach to the dynamics of large-disturbance 
maneuvers with particular emphasis on the nonlineari- 
ties associated with high rates of roll. We have tried 
to show what factors are important; (1) by comparing 
calculated responses, (2) by simplified equations and 
techniques, and (3) by study of simple models such as 
the force-free body. 

The problem is important and complex. 
that the designer is faced with now and will be faced 
For detailed calculations the com- 


It is one 


with in the future. 
plete equations of motion must be used and practical 
solution requires large digital or analog computers. 
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Present results are based on five degrees of freedom. 
However, even the five degree of freedom approach is 
open to question. It is obvious that maneuvers per 
formed at slightly greater than sonic speeds will result 
in deceleration into or through the transonic regime. 
The aerodynamic nonlinearities so introduced may have 
very large effects in rolling maneuvers and for such cases 
six degrees of freedom may be required. Other large 
disturbance motions are prevalent in aircraft which may 
be just as important as rolling namely, large pitching 


velocities and large forward accelerations. Studies 
have been made which show these too have destabiliz- 
ing effects similar to those of rolling. Thus, there 
is much work still to be done in the study of the dy 
namics of large-disturbance motions from the theoretical 
standpoint. 


There are a host of related problems that need solu 
tion. Aerodynamic data for use in large-disturbance 
calculations is meager, particularly in terms of non 
linearity with sideslip and rolling velocity. Schemes 
for automatically controlling the airplane in nonlinear 
maneuvers are needed, as are methods for design and 
synthesis of such systems. Design criteria are needed, 
not only in terms of the airplane control and structural 
limitations, but also in terms of acceptable rate and 
acceleration limits on the pilot. Even in the missile 
and ballistics field, where rolling has been a recognized 
problem for several years, there are many unanswered 
questions. 

However, the complexity of this problem in no way 
obviates the necessity for an elemental approach or a 
need for understanding its mechanics. Simple methods 
of approach are required not only for preliminary de 
sign studies but also for directing analog or digital cal- 
culations used in detailed design and analysis. Without 
such a guide the study of so complex a nonlinear system 


will become a hopeless hit or miss process. 

Various simplified approaches have been presented 
in this paper which show the effects of rolling velocity. 
These represent only a first attack of the problem. 
The overall domain of large-disturbance motions re- 
mains still primarily uncharted and the real need is for 
a clear understanding of these problems and simple 
methods of treating them which retain only essentials. 


APPENDIX 1 


EQUATIONS FOR TRANSFORMING FROM 
STABILITY TO Bopy AXES 


(1) Components 


xX = ¢ co a=- Zz sin @ 
Y= J, 
i= '¢ sin a + Z, COS @ 


The moments (L,, \/,, and .V,), the moment coefficients 


(Ci, Cms, Cns), the angular velocities (p,, g;, 7s), ete., 


all transform in similar manner. 
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? ‘electe aky ive _— . ° » : ‘ 
(2) Selected Derivatives Ihe transformation of derivatives where both quant; 
C,, = C,,,,c0ea — C,,,, ama ties are affected by the transformation requires a doub} 
. ‘ , ‘ operation. 
ct, = Cyy, 3 Sin avr ty , COS a 
Cin = Cyy,, COS a — (C; + Cny,s) X REFERENCES 
Sin a cos a + a sin” a 
. 4 os . . Rhoads, D. W., and O’Hara, J. C., Rolling Pullou tudy 
C, = C,., cos’ a + (¢ lp = Ge») om Part I, Development of Large Disturbance Theory for ¢ oe 
sin a cosa — Cus . sin? a Airplane Motions and Tail Loads ina Rolling Pullout Maney 
C. oe €.. AFTR-6701 Part I, Cornell Aeronautical Laboratory Repor 
é gree ' , TB-541-F-1, April, 1954 
Ce = &, cos YP (C.. . — Ce.) ‘ ; 
. iid ” aidnat ee x ‘ Te 2 Rhoads, D. W., Rolling Pullout Study Part Il, Flight 7 
sin a cos a — C;z,,, sin” a Veasurement of Vertical and Horizontal Tail Loads During a Rug 
Cc. = C,,,cos'a+(C;,.+ Cn ae der Locked Rolling Pullout Maneuver, AFTR-6701 Part II, Cor 
sin acosa+ C ie sin” a Continued on page 532 
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APPENDIX 2 


SUMMARY OF CONSTANTS IN EQUATIONS OF MOTION FOR THE SWEPT-WING FIGHTER 
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Optimum Nose Shapes for Missiles in the 
Superaerodynamic Region 


W. J. CARTER* 


lniversity of Texas 


SUMMARY 


The chanies of the kinetic theory 


on the nose of a missile moving in the super 


mic region of the atmosphere rhree separate cases are 


nsidered—ideal specular reflection, specular-type reflection 


lightly rough surface, and surface absorption followed by 


emission of the striking molecules. The calculus of varia 


en NOSE 


iployed to obtain the differential equation of the 
of the 


which minimizes the drag force for each three cases 


The resulting differential equations are then solved by a numerical 

cedure rhe drag coefficients for the optimum nose shapes are 
ikewise determined and these are compared with the drag co 
ficients given by other nose shapes. It is further shown that 
e drag coefficients arising when specular-type reflections occur 

significantly dependent on the nose shape When surface 
sorpti followed by random emission occurs, the drag co 
ficient is not strongly dependent on either the missile nose shape 
r the fineness ratio of the nose 


1) INTRODUCTION 


. VERY HIGH ALTITUDES, the mean free path of the 
A air molecules will be of the same order of magni- 
tude as the dimension of a missile moving through that 
region. When the mean free path is large, the drag force 
on the nose of a projectile may be computed, assuming 
that there is little mutual interference between the air 
molecules as they strike the nose and then rebound. 

In the classic formulation of the kinetic theory of 
gases, the molecules are imagined as being smooth, per 
fectly elastic spheres which will rebound from a smooth 
This 


The con- 


surface at the same angle at which they strike. 
type of behavior is called specular reflection. 
cept is obviously an overidealization of the matter in 
general. If all rebounds were of specular type, no vis 
cous type forces could be generated at a surface. Even 
though the classic picture is an oversimplification, it 
does account very well for many things. Experiments 
have shown that an appreciable percentage of specular 
type reflections can occur when a stream of gas molecules 
strike a smooth surface, especially when the impacts 
ire at an acute angle. 

For specular reflection, it is assumed that the dura- 
tion of impact is very short and that there is no oppor 
tunity for the striking molecules to achieve thermal 
equilibrium with the vibrating molecules on the surface 
of the solid. 
the rebound velocity is the same as that of the striking 


Under these conditions the magnitude of 


velocity. In such a case, the drag force may be ascribed 
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of gases is employed to 


to the momentum change of the rebounding molecules 


A slightly more realistic picture may be obtained by 
combining the concept of the smooth elastic sphere for 
the molecule with that of a surface having appreciable 
Fig. | 
The molecules rebounding from this surface are 


roughness. shows such an idealized rough sur 
face. 
considered to be individually of specular type and 
scattering occurs as a result of the surface irregularity 
Likewise, it may be seen that more molecules will strike 
the forward surfaces of the irregularities than will strike 
his action causes the average angle 


than that 


the rear surfaces. 
of reflection to be somewhat greater which 
would occur with specular reflection from the mean sur 
face. Some experimental evidence indicates that this 
sort of behavior occurs. 

Manv of the gas molecules which strike a suriace will 
be absorbed on the surface for a brief time interval and 
will then be emitted. The dwell time allows the gas 
the achieve 

the The 
velocity corresponding to the kinetic temperature, will 
Loeb 


cussion of this and other types of gas-surface phenomena 


molecule, on average, to thermal equi 


librium with surface. average emission 


have a random direction. gives a thorough dis 


2) OptimuM NOSE SHAPE FOR A MISSILE ASSUMING 
SPECULAR MOLECULAR REFLECTION FROM A SMOOTH 
SURFACI 


The curve of Fig. 2, when rotated about the x-axis, 


will form the nose of a missile. If the air molecules 
which strike the nose undergo specular reflection, the 
drag force on the missile may be deduced in a simple 
the 


occurs as the molecules rebound from the nose 


way by considering momentum change which 
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Nose shape of missile. 


The Eulerian point of view is taken here-—the missile 
is assumed to be at rest and the air molecules are con- 
sidered to be streaming past the missile nose in a direc- 
tion parallel to the axis of the nose. The magnitude of 
the streaming velocity is v. The air molecules also have 
thermal velocities corresponding to their mean kinetic 
temperature. It is usually assumed that the thermal 
velocities are random as to direction and 
Maxwellian distribution of magnitude. 

One of the significant parameters of free molecule 
flow is the ‘free molecule Mach Number,” 1/7, which 


have a 


has the definition 
Me = OE 


where c is the most probable thermal velocity of the 
molecules. Ashley* has shown that the thermal veloci- 
ties may be disregarded in determining the force ex- 
erted by the striking molecules, provided that v > c. 
For missile flight at very great altitudes this is the likely 
case. In the analysis which follows, it is assumed that 
M., is large enough so that only the streaming velocity 
need be considered. 

Referring to Fig. 2, it may be seen that the mass flux 
against the differential area 2aydy is 

dm = 2rpv ydy = 2rpv yy'dx 

where p is the mass density of the air and y’ = dy/dx. 


For specular reflection the rebound angle @ is the same 


as the angle of incidence. The total direction change 


relative to the axis of the nose is 28. 
The differential drag force due to the momentum 
change of the mass flux on the differential area is 


dD = 2mpv?(1 — cos 20)yy'dx 
from which the drag force may be obtained by integra- 
tion, 


D= apo? | (1 — cos 26)yy'dx (1) 
0 


It follows that 
1 — cos 20 = 2sin? 0 = 2y’2/(1 + y’?) 


thus, Eq. (1) may be written 
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D = 4n2pv" [ [yy’3/(1 + y’*) |dx 9 
JZ V0 

For the shape of the nose to be such that the drag 
force is minimum it is necessary that the integrand oj 


Eq. (2) satisfy the Euler-Lagrange differential equation 


(OF Oy) — (d/dx)(OF oy’ 0 3 


where F = yy’/(1 + y” A further requirement js 


that y and y’ be continuous and bounded in the interval] 


of integration. It may be shown that 


OF /dy = y8/(1 + y” 
OF /dy’ = [Byy’2/(1 + y””)] — [2yv’4/(1 + y’ 
(d dx)(OF Oy’) [6yy’y"/ (1 + y’”*)] + 
3y’3/(1 + y”)] — [l4yy’3y"/(1 + y’?)2] - 





[2y’®/(1 + y’*)?] + [Syyy?/(1 + y” 


Substitution in Eq. (3) and simplification yields the dif 
ferential equation of the minimizing curve 


5 fo 


y) [vw + 1)/(3 — y’?)] | 


y —(y” 
There is little hope of finding an analytic solution t 
Eq. (4) in terms of known functions, but a numerical 
solution may be readily determined. It may be seer 
that the differential equation has a singularity at y’ 

V3. For y’ < 7/3, the nose curve is concave down, as 
indicated in Fig. 3. For y’ > +/3 the curve of y is con 
cave up. It may be shown that y does not approach 
zero for any finite x. 
prefers to have a hole through the center. 


This means that the missile nos 
This is oi 
small practical concern since the point of curvature 
change occurs very near the nose centerline for nose 
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Comparison of nose shapes for //d = 2. 
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OPTIMUM 


NOSE 


The tip of 
1 to the x-axis may be 


shapes having a reasonable fineness ratio. 
the missile nose from say y’ 
formed by some suitable curve. In this case a circle 
are is taken to form the nose tip. 

The procedure for obtaining the numerical solution 
of Eq. (4 


mv convenient y; and y;’ as starting values (in this case, 
from Eq. 


may be outlined as follows by (a) choosing 
10 and y,’ 1.000); (b) computing y,” 
1): (c) choosing some appropriate interval Ax and 


computing 


Ye =V + yi Ax + y"(Ax?/2 


yo’ = yy’ + "Ax 


ind (d) repeating the process of (b) and (c) until the 
curve has been adequately mapped. 

Table 1 gives a few of the values of y, y’, and //d for 
the minimizing curve. A graphic comparison of the 
curve with a cone and a parabola having its vertex at 
the junction of nose and body is shown in Fig. 4. 

It is of interest to compare the drag coefficients at 
various values of / d for the three nose shapes shown in 
Fig. 4. The drag coefficient will be defined by the re- 
lation 


Cp 2D/A pv" (5) 


Eq. (2) gives D in terms of the nose shape. When this 


is substituted in Eq. (5) the result is 


Cp (32 /d*) [ [yy’3/(1 + y’*) Jdx (6) 
7) 

Eq. (6) has also been solved by a numerical procedure 
in the case of the optimum curve and by analytic 
methods for the cone and parabolic shape. These re- 
sults are shown in Fig. 5. It may be seen that the op- 
timum nose shape gives a somewhat lower drag co- 
efficient for all values of //d. 

A MISSILE ASSUMING 


MOLECULES 


3) OptimuM NOSE SHAPE 
REFLECTION OF THE 
SLIGHTLY ROUGH SURFACE 


FOR 


SPECULAR FROM A 


The Introduction has discussed the manner in which 
a molecule would rebound from a nonideal surface. It 
is further suggested that the smoothest of realizable 
surfaces have irregularities that are of the same order 
Most of the 
air molecules rebounding from the rough surface would 


of magnitude as a molecular dimension. 


undergo a change of direction greater than that con- 
sidered in Section (2). 

Some idea may ke obtained of the effect of surface 
roughness on the optimum nose shape by assuming 
that the molecules rebound from the surface as shown 


TABLE | 
\ y y’ 1/d x y y’ l/d 
+.14 10.00 1.000 0.207 500 152.7 0.208 1.64 
o 10.52. 0.949 (0.288 1,000 243.7 €.165 2.05 
13 16.60 0.645 0.392 3,000 517.7 0.119 2.90 
eV) 34.98 0.415 0.715 6,000 838.5 0.097 3.58 
100 53.50 0.335 0.985 15,000 1585 0.072 4.73 
720) 88.70 0.262 1.24 35,000 2836 0.056 6.17 
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in Fig. 6. The average rebound angle is assumed to be 
larger than the angle of incidence by the amount e. 
In the same manner as for the case of the perfectly 


smooth surface, the differential drag force is found to be 
dD 2rpv"|1 — cos (26 + e€)lyy’dx 


from which 


D 2m pv" [1 — cos (26 + e) |yy'dx 7 
7 


It may be also shown that 
cos (24 + e) cos 24 cos € — sin 26 sin ¢ 


To expedite the subsequent work, ¢ will be assumed to 


be of such smallness that cos e l and sin e = e—thus, 
cos (26 + e€) cos 26 — ¢€ sin 26 
D 2m pv" (1 — cos 26 + € sin 20) vy'dx s 
7 O 
It may be shown that 
1 — cos 26 2y’?/(1 + y” 
sin 20 = 2y’/(1 + y’ 
D 4arpv? X 
l { [yy’3/(1 + y’2)] + fevy’2/(1 + v2) ]fdx (9 
J V 


As in Section (2), the differential equation of the 
minimizing curve is found by substitution of the 
integrand of Eq. (9) in the Euler-Lagrange relation, Eq 


(3). A tedious bit of algebra discloses that 


‘. ar (Ch + y'*)* — el — 9" ] 

y —-* a pe : —— = 
: y [((38 — y’2)(1 + y’”?)? + €(2 — 6y”)] 

It will be seen that Eq. (10) reduces to Eq. (4) when « 


0. Eq. (10) has been solved by the same type of 
numerical procedure used in obtaining the solution to 
Eq. (4). 10 and y’ 
1.000. A circle arc is used to form the end of the nose 


shape and ¢ was chosen to be 0.07 radian. 


The solution is started at y 
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The roughness of the surface has very little effect on 
the optimum nose shape for e = 0.07 radian. This be- 
comes obvious by comparing the two curves shown in 
Fig. 7. It should be noted that the y-axis must be con- 
siderably stretched to show much difference in the two 
curves. 

Although the surface roughness has little effect on 
the optimum nose shape—at least for the mechanism 
it does have a significant effect on the 

The drag coefficient is defined in the 


assumed here 
drag coefficient. 
same manner as before: 


Cy» = (82/d*) X 


| riyy’? (1 + vy?) | telyy? (+ y%)]f dx (11) 

Jo 

Eq. (11) has been solved by a numerical procedure for 

the optimizing curve defined by Eq. (10), and the result 

is shown in Fig. 8 together with the drag coefficients 

that would result if the rebound angle equaled the angle 

of incidence—that is, « = 0. 

(4) OptimuM NOSE SHAPE ASSUMING MOLECULES ARE 

REFLECTED FROM THE MISSILE NOSE 
SURFACE 


RANDOMLY 


Experimental evidence indicates that many of the 
gas molecules which strike a surface are absorbed on the 
surface for a brief time interval and are then emitted. 
The percentage of gas molecules which will undergo 
this action is dependent on the temperature of the sur- 
face and the chemical nature of the gas and the surface. 

It is assumed that the molecules which are absorbed 
on the surface and are later emitted will leave the surface 
in a random direction with an average magnitude ot 
velocity corresponding to the surface temperature. 
The average velocity of emission is assumed to be 67, 
where v is the missile velocity through the atmosphere. 
Fig. 9 illustrates that random emission implies that the 
probability density of emission is uniform in the unit 
hemisphere above the point at which the molecule strikes 
the surface. The percentage of molecules emitted in 
the solid angle forined by the rotation of d¢ around the 


normal to the surface is 
(27 cos ddd) 27 = cos odo 


The velocity component of these molecules in the 


direction of the normal to the surface is 
dv, = bv sin ¢ cos ddd = (1/2) bv sin 2¢d¢ 


The net velocity component in the normal direction for 
all of the molecules emitted above the point is given by 
the integra] 
V, = (1/2) bv | sin 2¢d¢6 = bv 2 (12) 
70 

The conclusion derived from the foregoing analysis is 
that the average velocity vector of the emitted mole- 
cules is normal to the surface and that it has a magni- 


tude 


cv = bv/2 
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TABLE 2 

x ¥ y’ l/d ’ y y 

4.14 10.00 1.000 0.207 300 S2.58 0.157 t2 

7 2.35 0.709 0.283 650 129.3 O.11¢ ys 
12 15.46 0.552 0.388 1.500 211.7 0.083 3 
30 23.59 0.385 0.636 3,000 319.8 0.064 47 
60 33.91 0.295 0.885 $4,000 380.4 0.057 & 
130 51.50 0.218 1.26 5,500 461.5 0 051 


The mass flux on the dy strip is given by the expr 
sion: 
] =— 9 sda ») a ae 
am = pvva\ =—Tpvvy ax 
from which is obtained the differential drag force 


dD 


= 2rpv(v + cv sin O)yv'dv 


The total drag force is given by the integral 


D = 2xpv’ | (1 


J0 


+ c sin 6)yyv'ds 


It may likewise be shown that 
sin @ = y’ ‘(1+ y” 


» 


Thus Eq. (13) may be written 


D = 2npv? [ yy’ + feyy’? (1 +) *] de (1 
0 


JW 0 


is independent of the nose shape, and thus this part 


The integral 


vi ly 


of the integral of Eq. (14) need not be considere 
in the minimization procedure. The expression to b 
minimized in this case is 

»\ 1/2 

habs oS 


D = 2mrpvc | [yy’? (1 + y’ 
0 


For a minimum drag the integrand of Eq. (15) mus 
satisfy the Euler-Lagrange differential equation whic! 
vields the differential equation of the minimizing curv 
In this case, 


” 


y" os 


(1 + vy”) /(2 -— vy” 


— (y”? y) 


Eq. (16) has a singularity at y’ = +/2 which need cau 
no concern here. As before, the solution is started : 
y = 10, vy’ = 1.000, and the extreme tip of the no 
shape is formed with a circle are. Eq. (16) has bee! 
solved by the numerical procedure described in Sectio! 
(2). Table 2 gives values of x, y, y’, and / d for ths 
curve. The shape of the curve is very nearly that of 
parabola which has its vertex at the noce. 

It is convenient to define a coefficient C, which | 
dependent on D [see Eq. (15) ] thus, 


. 
(16c d?) [yy’? 
0 


(1 + y’?)°"|dx 


Eq. (17) shows that the portion of the drag c ye flicien! 
which is affected by the nose shape is also dependent 0! 
c which is a function of the surface temperature of the 
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nose or more exactly a function of the difference of the 
kinetic temperatures corresponding to the missile 


velocity v and the surface temperature of the nose. 


For the purpose of comparing the drag of various nose 
shapes with random reflection, the coefficient may be 


defined as 
/ 
Cp c (16 d?*) | [yy’? (1 + y’?)} i dx (18) 
J0 


This coefficient has been computed by a numerical pro- 
cedure for the optimum nose shape defined by Eq. (16). 
The coefficient has also been determined for a parabolic 
nose shape, having the vertex at the missile nose tip, and 
for a cone. A comparison of the values for these shapes 
is shown in Fig. 10. The values for the parabola are 
very near those for the optimum curve since the opti- 
mum curve has very nearly a parabolic shape. 

The total drag coefficients do not differ very much 
when random reflection occurs. The total drag co- 
efficient may be defined by the relation 


Cp = 2DAv’ 
(16 d?) | yy’ + feyy’?/(1 + y)'*]} dv (19 
/7 
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This coefficient has been evaluated for c equal to unit 
for the optimum curve and the cone. The results ar, 


— 


shown in Fig. 11. These curves show clearly tha | 


neither nose shape nor the / d ratio have much effec 
on the total drag when random reflection occurs. 


(5) CONCLUSION 


The fundamental information now available is not 
complete enough to permit one to determine exactl 
the aerodynamic drag of objects moving in a rarefie 
gas. As noted in the Introduction, the drag force may 
well be dependent on factors other than the geometr 
of the body. This study then should be regarded as on 
which defines the limits between which the actual drag 


forces must lie. 
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SUMMARY 


The equations obtained by Chien for the nonlinear deflection 
ff shallow spherical shells under uniform external pressure ar¢ 
solved by means of power series expansions, following procedures 
introduced by Friedrichs and Stoker in their treatment of buck 
ling of circular plates. These equations depend upon two param- 
eters. One of these parameters is related to the external pres 
sure, while the other depends upon the dimensions of the shell 
The 
under boundary conditions corresponding to a fixed edge 

The solution, carried out numerically on the AEC UNIVAC 
it New York University, yields a complete description of the 


stresses and deflections as functions of the polar angle over a 


equations are solved for several ranges of the parameters 


of values of the loading parameter and the dimen 


Prediction of the upper buckling load is then 


wide range 


sional parameter. 


made by means of a numerical criterion based on the load vs 
deflection curve. For some cases, the postbuckling behavior 
is investigated. The results agree well with existing experi 


mental and theoretical studies and cover a wide range of cases 


not previously treated 


SYMBOLS 


D = flexural rigidity of shell 2Eh 
Young's modulus 
boundary functions, see Eq 
= shell thickness 
= 2/3(1 v? 
= number of terms in power series required for con 
vergence 
external pressure, positive inwards 
= loading parameter, see Eq. (12 
radial variable, distance to any point in shell 
radius of middle surface of shell 
) = surface of shell 
= tangential component of displacement 
‘ = tangential component of displacement on middle 
surface 
= volume of shell 
= radial component of displacement 
= radial component of displacement on middle surface 
a = slope of tangent plane to deformed surface with the 
horizontal plane 
= &@&/A 
= stress function, see Eq. (10 
= reduced stress function, see Eq. (12 
longitudinal strain 
circumferential strain 
€ = asmall number 
components of strain tensor 
Presented at the Structures Session, Twenty-Fifth Annual 
Meeting, IAS, New York, January 28-31, 1957 


The research reported in this paper was performed under 


Contract AT(30-1)-1480 with the United States Atomic Energy 


Commission. A partial summary of the results contained in 


this paper was presented at the Ninth International Congress of 


Applied Mechanics in Brussels, September, 1956. 
The authors wish to express their appreciation to Prof. J. J. 


Stoker for suggesting this inv estigation 


or 


= a small number 


C7] = polar angle 

4 = 6/A 

A = semiangle of shell opening 

\ = Kaplan and Fung geometric parameter 

L = Poisson's ratio 

II = potential energy 

f = geometric parameter, see Eq. (12 

oj, 77 = longitudinal stress; value on middle surface 
o3; ¢3° = circumferential stress; value on middle surface 
T = components of stress tensor 

¢ = azimuthal angle 


INTRODUCTION 


pres a thin hollow sphere of wall thickness 
2h and mean radius R loaded by a uniform external 
pressure p directed toward the center. It is observed 
experimentally that as the pressure increases from 
zero, the shell apparently deforms by contracting into 
This mode of deformation continues 
At 


this value the shell suddenly jumps, or buckles into a 


a smaller sphere. 
until some critical value of the pressure is reached. 


nonspherical buckled state consisting of a round, rota 
tionally symmetric, inward ‘“‘dimple’’ of small solid 
angle at some point on the surface of the shell, while 
the remainder of the shell remains apparently spherical. 

These experimentally observed facts are entirely at 
Char 


acteristically, the linear theory predicts a wave form 


variance with the linear theory of buckling.! 


extending over the whole surface, with the amplitude 
e.g., the shell could buckle inward or 
outward load. to 
trasted with the experimentally observed ‘“‘dimpling”’ 
the 
load is found, in some cases, to be as low as one quarter 


undetermined 


under the same This is be con- 


behavior. Furthermore, experimental buckling 
of that predicted by the linear buckling theory. 

In a similar type of phenomenon observed in the case 
of buckling of thin cylindrical shells under axial com 
pression attempts were made to explain the discrep- 
ancy between theory and experiment on the basis of 
initial imperfections from the cylindrical shape of the 
model, and deviation of the boundary conditions from 
those the Al- 


though important, these effects were not sufficient to 


assumed in theoretical calculations. 
account for the lack of agreement (see the discussion in 
reference 2). 

Von Karman and Tsien* recognized that this must 
be the result of a nonlinear behavior of the structure 
Based upon the minimization of a potential energy 
expression, they showed that for values of p below the 


buckling load predicted by the linear theory there exist 
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h= = thickness 


Fic. 1 Shell geometry 
other, nonadjacent, stable states of equilibrium. The 
greatest lower bound of such loads /p is called the 
“minimum” or “lower” buckling load px. If p exceeds 
Px there is a possibility that buckling will occur. 
Friedrichs*® pointed out some shortcomings and sug- 
gested improvements in the Karman-Tsien theory. 
He also treated the buckling of a thin sphere as a 
“boundary-layer? phenomenon; where the ‘boundary 
layer’’ is considered as occurring in the dimpled region 
of the sphere. Tsien,’ adopting some of the sugges 
tions of Friedrichs,* presented an improved version of 
the Karman-Tsien’ theory. Altshuler? attempted 
to determine the buckling load by solving the differ 
ential equations by asymptotic expansion techniques. 
Although none of these approaches are entirely satis- 
factory, they do reveal the essential nature of the phe 
nomena that are involved in the buckling process. 
Closely related to the problem of the buckling 
thin sphere is that of the buckling of a shallow* spher- 
Physically 


fa 


ical shell under uniform external pressure. 
as well as mathematically, both problems reveal a 
similar nonlinear behavior. The shallow spherical 
shell is of interest both as a structural entity in itself, 
and as representing that portion of a complete sphere 
in which the dimpling occurs. 

Solutions for the nonlinear deformations of shallow 
spherical shells under various external pressures and 
edge moments have been obtained by Chien and Hu'” 
and Feodosyev.'* In reference 6, Kaplan and Fung 
present a theoretical and expertinental treatment for 
the clamped shallow spherical shell under uniform 
external pressure. Their theoretical work is based 
upon a perturbation solution of the nonlinear equa- 
tions, where the nondimensional radial deflection at 
the center of the shell serves as the perturbation 


parameter. Unfortunately, solutions to the resulting 


By shallow it is implied that the semiangle of opening A 


of the spherical segment is a small angle (see Fig. 1 


‘ 


linear perturbation equations could be obtained on] 
for very shallow shells. However, their experiments g 


observed in reference 7 revealed some new and inter 


esting phenomena concerning the buckling of curve 


structures. 


derived by Reissner,'® extended and improved th 


perturbation technique, and was able to obtain soly 
tions for a greater range of shells than Kaplan an 
Fung. The agreement with experiments of the buck 
ling loads predicted from the perturbation  solutio; 


becomes poorer as the value of the geometric shel] 


A*R Vkh) increases. Furthermor 


parameter (p 


the perturbation solution, based upon a one parameter 


expansion, fails to reveal the nonmonotone dependence 
of buckling loads on p as indicated in reference 7 
Simons'? proposes solving the nonlinear equations"® by 
a power series technique similar to that presented i1 
this paper. His calculations, however, cover only 
small range of parameters. 

The physical problem treated in this paper is that 
a thin shallow spherical shell under the action of um 
form external pressure. 
clamped along the edge. Only rotationally symmetri 
deformations will be considered. The problem is con 
sidered as one in the nonlinear theory of elasticity wit! 


small strains and small, but finite, deflections. Non 


linear load vs. deformation characteristics of the shell 


are expected to have the form shown in Fig. 2. Here 


p represents the external pressure and w the radial de 
flection at some point—-i.e., the center. Point A rep 
resents the “upper” or “initial” buckling load. Point 
B represents the ‘final’ or ‘‘minimum’’ buckling load 

The present investigation has been successful in ob 
taining solutions including deflections and stresses fot 
a large range of shell geometries for pressures fron 
zero to those in the vicinity of point A of Fig. 2. Ur 
fortunately, only meager results have been obtaine 
so far for the “postbuckling”’ region —i.e., the portior 
of the curve ABC in Fig. 2. 

The technique of solution employed was that of power 


series expansions inspired by the success of Friedrichs 
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NONLINEAR DEFLECTIONS OF 


and Bromberg! in the problems of bend- 


and StLOKCT, 
ing buckling of circular plates. The use of this 
approach in the present problem, however, was only 


sossible due to the availability of a high-speed digital 


computer 


DERIVATION OF DIFFERENTIAL EQUATIONS 


The nonlinear equations for the deformation of 
spherical shells of small curvature, based upon an 
issumption of small strains and linear stress-strain 
law, were originally obtained by Chien* from static 
considerations. These equations can also be derived 
from a variational principle. 


the principle of minimum potential energy and seek 


In this case we employ 


Euler's equations for the functional, 


¥ 
a 
‘ 
> 
a 
> 
~™ 
> 
S 
ce 
& 
S) 


where II is the potential energy, 7” are the components 
of the stress tensor, €,; are the components of the strain 
tensor, w is the radial displacement, and / is the ex 
ternal pressure, positive inwards. See Fig. 1 for the 
shell geometry. 

The functional, Eq. (1), may be expressed in terms of 
the components of displacement through the following 


relations. 


Stress Strain 


From the assumption of rotational symmetry all other 


stress and strain components vanish. 


Strain Displacement 


Let r be the radial distance to any point in the shell 
and R the radius of the middle surface. The exact 
nonlinear strain displacement relations for rotationally 
syminetric deformations in a spherical coordinate sys- 


tem are then 


ucotd — W)/f| + 
1 2) [(u cot@ — w)/r]? (3b 


Where subscripts on « and w denote differentiation. 
For the case of small finite deflections, « and w (anal- 
agous to the Foppl-Karman nonlinear plate theory) 
and shallow shells (small 6, uj? is negligible compared 
to uj and w,*, and uw is negligible compared to wy), 


Eqs. (8a) and (3b) reduce to 


€j (l/r) (ug — w) + (1/2) (w6/r ta 
es = (1 r) (u 0 — w) (4b 
Eqs. (4a) and (4b) are the strain displacement rela- 


ions that will be employed in the variational expres- 
sion. It is not clear to what extent these equations are 
adequate for nonshallow shells without retaining addi- 
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tional terms from Eqs. (3a) and (3b); for example, see 
page 260 of reference 3. 

As is customary in thin shell theory, the following 
assumptions are made concerning the distribution of the 


displacements through the shell thickness: 


u u°—(C Ruz: «z Tp 5 
where ¢ R — rand u’ and w’ are the displacements 
of the middle surface ¢ 0. 

Substituting Eqs. (2), (4a), (4b), and (5) into Eq. 


(1), and performing the ¢ and ¢ integrations, the po 
tential energy II may be written as 


vA 
11 /27R? [Hh (1 — 1 €;°* + «3, 4 


A 
2veEz es") 6d@ + (FE: | x 


(h® R*) [wag + (1 6?) wz? + (2v O)wws;) 6d6 4 


7A 
bwhdé 6 


where ¢;’ and e;° are the strains on the middle surface 


Setting the first variation of Eq. (6) equal to zero, 
the following equilibrium equations result: 
os; = d d6) (80; a 


Sy \ ) ( 
— (2h/R) }(o;° + o3°)0 + (1 R) (d dO) (Ow 5 0; } t 


(D R*) j(d d6) 6(d dd) (1 6) (d db)bw;} — pd = 0 (7b 
where o;° and o; are stresses on the middle surface, and 
D PEA BL — vy 


Eqs. (fa) and (7b) are a set of two equations in the 


three unknowns, o;’, o;’, and w. It is also necessary 
to satisfy the compatibility equation for the middle 
surface. This equation can be obtained by eliminating 


u from Eqs. (4a) and (4b) forr = R. There results 


(d db) 60; — we) - FF + Ws 


(E’'R)6O wz + (E 2 


(Ws R -= WU (S 


In addition to yielding the differential equations 
(7a) and (7b), the variational problem leads to the 


following boundary conditions: 


[(2hR)o;'06u" + (6 'R2) X 
\(D R2) [—(d dé) (16) (d db)bws + vws] + 
2hoiw;t é6w + (D R4 (Oza + vwy;) dw, = QO (9 


The first equation of equilibrium, (7a), is identi 
cally satisfied if we set 

oj (E 2h) (7/6 (10a 

a! E. 2h) 7; 10b) 


Also, let us introduce the more convenient variables 


a (1 R dw d6) + 6 (11) 
a Aa 6 = A6 k 2/3(1 — 2 
7 k (Ah? R)y; p 2E(R h/R)*P: (12) 


p k)- PAR I 
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28 : 


Fic. 3. a, y; plane. 


Geometrically, the independent variable & defined 
in Eq. (11) is the angle that the tangent plane to the 
deformed surface makes with the horizontal plane, see 


Fig. 1. 


The parameter p, hereafter referred to as the ‘‘geo- 
metrical parameter,’’ is proportional to the ratio of the 
shell rise and the thickness. The shell rise is defined 
as the height of the center of the shell above the hori- 
zontal plane through the edge. p is also roughly pro- 
portional to the ratio of membrane and bending stiff- 
ness. The relation between p and Kaplan and Fung’s 
geometrical parameter ) is 


p = (2)—"/*)2 (12a) 


Substituting Eqs. (10a), (10b), (11), and (12) into 
Eqs. (7a), (7b), and (8), and integrating Eq. (7b) once 
(employing the fact that from rotational symmetry 
dw/d@ = 0 at 6 = O), the following differential equa- 


tions are obtained: 


La = play + Pé?) (13a) 
Ly = p(@? — a’) (13b) 

where L is the linear differential operator, 
L = 0(d/d6) (1/0) (d/d0)6 (14) 


The boundary conditions corresponding to the dif- 
ferential equations (13a) and (13b) are as follows. At 
the center, @ = 0, from rotational symmetry, 


u=0 


From Eqs. (2), (4a), (4b), (10a), and (10b), this is 
equivalent to 


a(O0) = 0; y(0) = O (15) 


SICAL SCIBNCES FUL Y, 


1957 
At the edge, 6 1, we assume that the plate is hel 
fixed and clamped. 
as simulating the boundary separating the 


This condition can also be cop. 
sidered 
region of dimpling from that of uniform contraction oj 


a complete sphere. Therefore, the conditions ar 


This is equivalent to 
a(l) l; [yo — vyleus 0) lf 


It should be observed that, if conditions (15) and (16 
are imposed, then the boundary conditions (9), resulting 
from the variational problem, are identically satisfied 

Thus, the complete formulation of the problem con 
sists of the differential equations (13a) and (13b) subject 
to the boundary conditions (15) and (16). Theg 
equations——-a pair of second-order coupled ordinan 
nonlinear differential equations 
the equations governing the nonlinear deflections oj 
'! except for the addition of extra 


are quite similar t 


circular plates,° 
terms resulting from the initial curvature. 
the physical phenomena that these equations describ 
are quite different. Mathematically, the buckling 
problem for the circular plate, or the bending problem 
differ in their parametric dependence from the spherical 
The plate problems depend upon one 


I I< ywever 


shell problem. 
parameter, either a reduced edge load (the ‘“‘thrust 
ratio,’’ see reference 10) or a reduced lateral load (see 
reference 11); while, on the other hand, the equations 
governing the shell problem depend upon two param 
eters, P and p. From this point of view, it would seem 
that to obtain a more reasonable perturbation solutio 
two-parameter expansio1 


for the shell problem a 


should be attempted. 


TECHNIQUE OF SOLUTION 


The technique of solution to be employed for the 
boundary value problem posed by Eqs. (13a), (13b 
(15), and (16) is to expand the dependent variables « 
and y in power series in @. This approach has beet 
used in nonlinear circular plate problems by Way 
Friedrichs and Stoker,'® and Bromberg.'! 

From Eqs. (13a) and (13b) and the boundary cond 
tions at 6 = 0, the power series for a and y must con 
tain odd powers only, so that 


a= > a,0""; y= > 7,6?" (17 
n=1 


By substitution of these power series into Eqs. (19a 
and (13b), the following recurrence relations for the 


coefficients a, and y, are obtained: 


a. = p/SX(ay+ P); ve = p/8(1 — an’); 
for n > 2 
p " 5 IS 
Qa, = a1Yn-i; , 
‘ 4n(n — 1) p> si 
n—1 
p t 
= — AiAyn 
4n(n — 1) p> 
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te is held Again, by insertion of the power series into Eq. (16), TABLE | 
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action oj 
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j condi: 


ist con 


for the 


the boundary conditions at 6 1 can be written as 


> a, 1; >} (Qn —1— )y, = 0 


— 
1 ] 


Since it is possible through the recurrence relations 


18) to express all a, and y, in terms of a; and y, the 


above conditions constitute a pair of simultaneous 


equations in a, and y;—-1.e., 


I a, Y1 » & a ics 0 
D 1 
(19 
G ay, Y¥1 Zz |2n _ L +% Yn 0 


l 


[he problem thus reduces, for given P and p, to solv 
ing the two simultaneous equations (19) for the un- 
knowns a; and y;. Once the root or roots of these 
equations are obtained all the a, and y, and, hence, all 
the stresses and displacements, are then determined 


see Appendix A). 


COMPUTATIONAL TECHNIQUE 


The computational problem is, for a given set of 
parameter values (/’, p) to determine the roots (a1, 7; 
satisfying Eqs. (19). For computing purposes, it is 
necessary to replace the exact Eqs. (19) by the approxt- 


mate conditions, 


Fla Y1 = > B a — ] | 
- (19’ 
Gla, ¥1 = z. |2n — (1 + V)\¥n| <1 


where 7 is a small number and NV is an integer deter- 
mined on the basis of a convergence criterion for the 
series. In the majority of the computations 7 was 
chosen as 1 X 10 

For a fixed P and p the procedure starts from a rea- 
, say (a1, yi1)°. Then from Eq. 


.V is never greater than 99. 


sonable guess for (ai, ¥; 


18) a, ® are computed 


and y, corresponding to (a1, y1) 
forn < N. As a convergence criterion for the series 
(17), it is preseribed that the sum of the absolute mag- 


nitudes of the last three computed a, and y, be less 


than a certain small number, «. That is, the com- 
putation of a, and y, in Eqs. (18) ceases for n N 
such that 
AN or QAN~1 — an? ' e/ 
(20) 
wir |trw-ti F Yn-2 < €\ 


In the majority of the computations « was set equal 
tol X 1074, 

When the a, and y, satisfy Ecs. (20), they are in- 
serted into the expressions for F(a, y:) and G(a, 1). 
li F and G satisfy Eqs. (19’), the stress and displace- 
ment distribution throughout the shell are computed 
from the relations listed in Appendix A. In general, 
however, F and G do not satisfy Eqs. (19’), and it is 
hecessary to obtain an improved approximation to 
a, 71). This is done systematically by means of the 


The Number JN, of Coefficients Required for Convergence 


P 
p 0.2 OS 1.0 1.4 
12.23 1] 13 14 19 
l7 .22 13 17 19 23 
32.0 16 26 31 
52.2 20) 36 1} 
64.0 22 42 a4 


Newton-Raphson iteration process. The partial de 
rivatives of F and G in the Newton-Raphson process 
are approximated by forward difference quotients. 
This iteration process is repeated until Eqs. (19’) are 
satisfied. The displacements and stresses are then 
computed from the equations of Appendix A. The 
number of iterations required for convergence varies 
with the (?, p) selection and the initial guess (aj, y;)' 
For most of the cases encountered, two to six Newton 
Raphson iterations were necessary. The number of 
coefficients, NV, that were required for convergence 
varied with the (?, p) selection; Table 1 gives some 
The 


process, 


representative values. computation of coeffi 
cients, the iteration 
stresses and displacements are all performed on the 
AEC UNIVAC at New York University.* 


The method for calculating with fixed ? and p has 


and the calculation of 


thus been described. The computations proceed, for 
fixed p, by obtaining a sequence of roots (a), y) for a 
corresponding sequence of increasing values of P?. A 
first guess for (ai, y:) is determined by selecting a low 
value for P, usually P = 0.1. 
P, the deformed state, and hence the first (a), 7: 
The Newton-Raph 


For such low values of 


guess, should be nearly spherical. 


Weissner and Messrs 


uid in coding and pro 


* The authors are indebted to Miss M 
A. Paster and S. B. Weinberg for their 


graming this problem 
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Fic. 5a Reduced radial deflection, p 11.312, Region I 

son process then rapidly yields the root (ai, y:) to the 
required accuracy. The initial guess of the roots for 
successively higher values of P is obtained by extra- 
polation of the preceding roots in the sequence. This 
is done with the aid of a plot of a; vs. y; as shown in 
Fig. 3. This procedure is terminated for either of two 
reasons—the Newton-Raphson iteration process fails 
to converge to a root in a large number of iterations 
(usually 10), or the coefficients fail to converge accord- 
As an example of the first reason for 
termination, with p = 32 and P 
However, with P? Be: 


ing to Eqs. (20). 
1.25 the procedure 
converged in 3 iterations. 
there was no indication of convergence in 10 Newton- 
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Fic. 5c Reduced longitudinal total stress—top surface 


Raphson iterations. An illustration of the second 


reason for termination 1s for p 26.1 and P 1.55 


Here convergence occurred in 2 iterations and 40 terms 


But, for P 1.6, 99 terns were 


computed with Eqs. (20) not being satisfied. 


of the power series. 
In fact 


the terms were steadily growing in absolute magnitude 


The criterion for the initial buckling that is used in 


this investigation is based on this terminating value 
of F. 


(the load at which dP dw = O) for buckling seems im 


Direct application of the “‘classical’’ criteriot 
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practicable with the present technique of solution 
Instead, the numerical equivalent of the buckling 
phenomena is used as the criterion for determining the 
buckling load—.e., 
ess to produce a neighboring solution corresponding to 


the failure of the computing proc 
small increment in the load. If 7; represents the 
last value of P tor which convergence occurs, and 7? 
presents the next value of P? attempted, then the 
initial buckling load, P 
terval Py. Per SS Ps. 


is assumed to he in the in 


DISCUSSION OF RESULTS 


Before proceeding to a discussion of the results of 
our calculations, a brief review of the experimental re 
sults of Kaplan and Fung and the interpretation given 
in reference 7 will be presented 

Kaplan and Fung constructed a series of shallow 


These 


spherical domes of varying radii and thickness. 
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Fic. 7b. Reduced longitudinal stre 32, Region II 


domes were placed in a clamping jig and loaded under 
uniform external pressure. In addition to measuring 
the initial and final buckling pressures, they measured 
the vertical deflections at various points for several 
different values of the pressure. 

For low values of the geometric parameter p [ort 
equivalently Kaplan and Fung’s \ as related by Eq 
l2a)| the experiments showed that the form of the 
deflections are such that the maximum deflection is at 
the center of the shell, and decreases steadily toward 
the edge. This will be called mode of deformation 
I. As p increases, it was observed that the pe ik de 
flection at the center gradually flattened out, until, at 
p equal to approximately 21, the maximum deflection 


occurs at some point between the center and edge 
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Reduced longitudinal total stress 
p = 382, Region IT. 


Fic, 7d. 


so that there is also a local minimum point at the center. 
With a 
further increase in p the maximum point moved toward 
the edge until approximately p = 55. Then the 
shell deformed in a third mode. This mode has two 
local maximum points, one at the center and the other 


This will be called mode of deformation II. 


between the center and the edge (see Fig. 16 of refer- 
ence 6). 

When the experimental buckling loads are plotted 
against p rather than X, it is found’ that there is an 
interesting peaking behavior in these loads as a func- 
tion of p—e.g., see Fig. 4. In Fig. 4 the crosses and 
circles indicate the experimental values of Kaplan and 
Fung, while the solid line indicates the apparent trend 
of the experimental data. It is significant to note that 
these peaks in the buckling loads occur at those values 
of p at which a change of mode of deformation had been 
observed. 

In proceeding with our numerical solutions on the 
UNIVAC we have in a real sense conducted ‘‘experi- 
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Fic. 8. Reduced radial deflection, p = 176, Region III. 


bottom surface, 


ments’? on shallow spherical shells. For each give 


value of p—i.e., a particular shell configuration 
deflection and stress 
increasing values of load P. As in a physical exper; 
ment, the initial buckling load is considered to ly 
reached, when for a small increment in P no neigh 


boring solution to the previously attained solutioy 


can be found—i.e., the iteration process suddenly fai 
to converge. 

For low values of p, p < 23, the numerical results 
confirm Kaplan and Fung’s experiments. In Figs 
5a, 5b, 5c, and 5d, are plotted some values of the de. 
The plots are 
for increasing load P at a constant value of p. In this 
11.312 was selected as representative of the 


flections and stresses as a function of @. 


case p 
mode of deformation of type, or region I. 
the radial deflection w can be seen to peak at the center 
and fall off toward the edge. The dotted lines indicat 
ing the results from the perturbation solution of Kap- 


In Fig. 5a 


lan and Fung.’ It can also be noted from Figs. 5b, 5c, 
and 5d that the stresses exhibit a similar behavior. 
For increasing values of p, the graphs of these func- 
tions tend to flatten out in the central region of the 
shell (9 = 0). 


See Fig. 6, where only the results for the radial de- 


23.5. 


This process continues until p 
flection are presented. For this value of p and for low 
values of load, say P < 0.4, the shell deforms in mode 
I with maximum deflection at the center. However, 
for higher load the mode of deformation gradually 
changes into type II, with the maximum deflection 
located somewhere between the center and the edge 
(in this case, 6 ~ 0.3). It is to be noted from Fig. 4 
that the shell p = 23.5 falls in the peaking region of the 
buckling loads. Thus, the mode of deformation for 
constant p may change with increasing load. This 
clearly illustrates the two parameter dependence of the 


character of the solutions. 











Reduced radial deflection, at P = 1.0 
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7a, 7b, 7c, and 7d are shown plots of de- 


32, of mode 


fections and stresses for a typical shell, p 
Notice 


Il of deformation for increasing values of load. 
that the longitudinal membrane stress of Fig. 7b also 
has a maximum located between the center and edge 
5b). The the total 


(see Figs. 7c shows a 


compare to Fig. behavior of 
stresses for 
marked difference from that of mode I (see Figs. 5c 


For mode I, the total stresses fall off in a 


mode II and 7d 


and 5d 
smooth manner from the center to the edge, while for 
0.6 and then 


This is 


mode II there is a local minimum near 6 
the stresses change rapidly toward the edge. 
reminiscent of the boundary-layer phenomena found 
by Friedrichs and Stoker!® for the buckling of circular 
plates. 

A further increase in p results in a flattening of the 
deflection and stress distributions in the central por- 
tion of the shell, and a movement of the maximum 
point toward the edge. Mode II persists until p A 52, 
when a transition case similar to p 
However, this time the transition is from mode II to 
mode III. In the radial deflection plot, the mode III 
form is barely visible, but it can clearly be observed 
It is not until large p, say 


23.5 is detected. 


in the total stress graphs. 
176, that there is a sharp appearance of mode ITI 
Again 


p 
in the radial deflection graph (see Fig. S$). 
from Fig. 4 it can be observed that the transition of 
modes occurs at the peaking of the buckling loads. 

In Figs. 9a, 9b, 9c, and 9d the deflections and stresses 
are again plotted against @, but this time for constant 
P and increasing p. In Fig. 9a the motion of the maxi- 
mum deflection point toward the edge with increasing 
This is again reminiscent of 
boundary-layer the for the 
longitudinal membrane stress of Fig. 9b, and in par- 


p is Clearly indicated. 
phenomena. In plots 
ticular the total stresses of Figs. 9c and 9d, the change 
in mode shape is clearly illustrated. 

The load vs. deflection at the center characteristics 


for various values of p are plotted in Fig. 10 and com- 
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pared with the theoretical and experimental results 
obtained by Kaplan and Fung. 


the agreement with the experiments is only fair, while 


For low values of p 


for increasing p the agreement between the physical 
and numerical experiments becomes uniformly better. 
This is probably caused by the fact that for shallow 
shells (small p) the effect of initial imperfections is 
more pronounced, and that the shallower the shell the 
more difficult it is to measure the radius properly. 


Initial Buckling Loads 


The results of the numerical computation of 7P,, as 
well as the experimental results of Kaplan and Fung are 
presented in Fig. 4. The intervals, ?)/?2, in which the 
P., lie are indicated by short vertical lines. 
observed that the computed P,, follow the same gen 
eral trend as the experimental data, indicating a peaking 
The agreement between 


It can be 


behavior as a function of p. 
experiment and theory is least satisfactory in the re- 
gion of peaking [IS < p < 25, 50 < p< 60]. In this 
region a small error in the experimental determination 
of p may lead to a large error in P., for a particular 
shell. It can also be observed from Fig. 4 that in some 
cases the numerical results predict greater ?,, than the 
experiments, while in other cases it 1s just the opposite. 
There is no apparent trend in this respect. 


POSTBUCKLING BEHAVIOR 


The postbuckling region is defined as the region 
ABC of Fig. 2—1.e., the behavior of the shell after 
F F The primary purpose of investigating this 
region is to determine the lowest value of P-——(P 
P,)—1tor which there exists non-unique solutions to the 
boundary value problem described by Eqs. (15a), (15b), 
(15),and (16). The procedure previously outlined does 
not enable us to continuously obtain solutions beyond 
point A. However, starting from a solution in the 
neighborhood of point C we may apply the previous 
procedure with a decreasing sequence of P? values to 
obtain the point B. 

To obtain a solution in the neighborhood of point 
C we select a sufficiently /arge value of P (say P = 2.2). 
For a sufficiently sma// value of p it is possible to ob- 


TICAL 
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tain a solution, using the previous procedure, for thj 
value of P, since, for a ‘‘flat’’ shell, the load vs. defo, 
mation characteristics must be monotone—i.e., ther 
will be no buckling. Then, starting with the solutio 
ior this value of p, we obtain a sequence of solution 
for an increasing sequence of p values and fixed P (se 
Fig. 11; p 0.8, 4.8, 6.4, and 7.24 This procedur 
is analogous to the original procedure where the role: 
of P and p are interchanged. We note how, for j 
creasing p values, the characteristics depart more an 
more from the initial linear relation and approach th 
buckling type of curve. 

A solution is thus obtained for P 
11.312. This load is higher than that previously use 
and corresponds to much larger deflections. Keeping 


p fixed, a decreasing sequence of P values 1s selected i1 


2.2 and p 


an attempt to locate point B. Unfortunately, con 
vergence of the computations becomes increasingh 
poor as P decreases. Solutions are obtained only for 
the range 1.8 < P < 2.2, see Fig. 11. 

11.312 is also shown in Fig 


The initial por 
tion of the curve for p 
11. Fig. 12 represents the radial deflection over th 
shell for one value of ? in the postbuckling region 


CONCLUSIONS 


Numerical solutions of the nonlinear equations 
governing the deflections of clamped shallow spherical 
shells under a uniform external pressure have been ob 
tained by power series technique. The computations 
were conducted on the AEC UNIVAC at N.Y.U 


These solutions yield stress and deflection distributions 
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NONLINEAR DEFL 
over the shell as well as estimates of the initial buckling 
joads for a wide range of values of the shell parameters. 
The two parameter dependence of the solutions and 
the various modes of deformation have been clearly 
shown. Only a limited number of solutions for the 
postbuckling region have been obtained. It is an- 
ticipated that some improvements in the numerical 
method may yield additional results in this region. 


APPENDIX 
Once the roots (a), y;) are obtained for a given set of 
parameter values (7, p), the stresses and deflections 


may be computed from the following relations: 


\ 
a(@) z. a,,07" A-1) 
n=l 
wi) phVk (1/2),;1-— >: (a, n) + 

n=l 

\ 
> (a, ne" — 62} (A-2) 

=I 


pivVk (1 2)/1- 5 (a, n | A-3) 


w(O 
l 
Membrane stresses: @ o/F, 
\ 
2(k) "(Rhys = Do y,6?" A-4 
1=1 
\ 
2(Rk (Rh)é” = >> (2n — 1)7,62"-2 (A-5 


] 


lotal stresses: at upper and lower surfaces, 


\ 
R/h)&o\ = Vk 2)> > ¥,8 ht ms x 
wm 
| > n= 1 eS Wee” * = 1 se | (A-6 
\ 
Rhjé (V/ k/2) > (2n — 1)y,6 t 
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dimensionless variables into the heat-transfer equations 
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transition based on the distance to the cone tip divided by the 


plate Reynolds Number for transition based on the distance to 


the leading edge is not, in general, equal to three. The indication 
is that the ratio of transition Reynolds Numbers will be near three 
only when transition occurs near the minimum critical Reynolds 
Number When transition occurs at a much larger Reynolds 
Number, the indication is that the ratio will be near unity. Be 
cause these remarks apply only when transition is caused by the 
amplification of very small disturbances, they are not applicable 
when local surface imperfections or distributed roughness cause 


transition by generating large disturbances 
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T A PREVIOUS NOTE in the Readers’ Forum,' Sears function * 
S(w) was tabulated for » up to 15. As was pointed out in 


reference 1, this function (first used by Sears?) can be identified as 


S(w) = (2/mw)}1/[Jo(w) — Vilw)] — t[Silw) + Yo(w)] 5 | 


* The same function was denoted by F(w) in reference 2 and ¢(f) in ref 


erence 3 
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Fic. 1. Function S(w) 
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TABLE 1 
Function S(w 


wd Re (S) Im(S) “w Re(S) Im(S) 
0 1.00000 0 26 +-07793 +. 00678 
1 +. 36865 +. 12594 27 +.03573 +. 06795 
2 +.08157 +.26797 28 -.03719 +. 06559 
3 -.14479 +-17727 29 -.07403 +. 00341 
4 -.19802 +-02067 30 ~.04266 - .05904 
5 -.08116 -.15863 31 +.02620 - .06671 
6 +.08127 -.14098 32 +.06915 - .01380 
7 +-.15050 -.00765 33 +.04883 +- 05004 
8 +.08235 +.11443 34 -.01575 +- 06657 
9 -.04863 +. 12373 35 -.06360 + .0224) 
10 -.12365 +.02477 36 -.05246 - .04081 
ll - .08367 -.08635 37 +- 00587 - .06534 
12 +.02618 -.11212 38 +-.05735 - .03000 
13 +.10421 -.03712 39 +.05550 +.03162 
14 +.08439 +.06512 40 +.00335 +.06301 
15 -.00882 +.10261 41 -.05061 +. 03639 
16 - .08820 +.04654 42 -.05647 - .02345 
17 -.08420 -.04760 43 -.01175 -.05970 
18 -.00539 -.09390 44 +.04355 -.04170 
19 +.07404 -.05380 45 +.09785 + .01384 
20 +.08312 +.03240 46 +.01937 +.05554 
21 +.01728 +.08530 47 -.03589 +.04582 
22 -.06101 +.05930 48 -.05715 ~ .00536 
23 -.08101 ~.01887 49 -.02627 -.05056 
24 -.02730 -.07670 50 +.02815 - .04888 
25 +.04875 -.06311 


Recently, it was found that the values of S(w) up to e 
50 or more (see below) were required 
Thomas Water Tunnel made such a calculation 
are shown in Table 1. Since the accuracy of the values give 
in reference 1 was found to be insufficient, S(w) for w up to 15 wa 


recalculated. In so doing an error in the previous reference (for 


w = 10) was discovered. This was corrected in the present table 
For values of w larger than 50, S(w) may be obtained by t 

following approximate relations 

Jo(w) > V 2/mw cos [w — (2/4 

Vo(w) > VY 2/2w sin [w — (2/4 

Ji(w) > VY 2/7w cos [w — (32/4 

Vio) ~¥ 2/rw sin lw — (37/4 
Using these relations, it follows that, for large values of w 
S(w) > (1/2 Y zw) [(sin w + cos w) + 7 (sin w — cos w)| = 

1 — 1)/2V role ~ ia 


S(w)| > 1/*V 22w 


Eq. (3) was first suggested in reference 3 (which was concerned 
with the power spectrum and therefore required only the magni 
tude but not the complete calculation of this function \ com 
parison using the exact formula (1) and the approximate formula 
2) shows that the discrepancy is negligible for w > 50 

Eq. (2) not only facilitates the calculation of S(w) for larg 
values of w but also points out the now obvious fact that, for large 
values of w, S(w) rotates in the same manner as w—i.e., S(w) makes 
one revolution in the vector diagram (see Fig. 1) when w increases 
by 27. (A similar diagram for 0 < w < 10 was given in reference 
2 Furthermore, Eq. (3) shows that the magnitude of S(w) de 
creases slowly with increasing w. Physically, this means that the 
response of an airfoil to a sinusoidal gust drops off only gradualls 
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READI 
ing frequency of the unsteadiness. This is the reason 
ip tow 50 or more may be actually necessary in 

ation of a practical problem (such as the problem 
f pressure perturbation due to sinusoidal upstream velocity 


It is hoped that Table 1 may prove helpful to 


jistul ‘ 
the estigators in similar problems 
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The Pressure Rise at Shock-Induced Turbulent 
Boundary-Layer Separation in Three- 
Dimensional Supersonic Flow* 
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April 1, 1957 


University 


SUMMARY 
Experiments are reported which indicate that the peak-pressure ratio for 


ersonic turbulent free separation depends only on the Mach Number 


onent normal to the ne of separation 


EXPERIMENT 


K* ERIMENTS Were undertaken to study the three-dimensional, 
shoek-induced turbulent boundary-layer separation process 
it supersonic speeds. Some results obtained regarding the peak 
pressure ratio attained in fully developed or “‘free’’ separation 
are sufficiently important to warrant immediate publication as 
this mterim note 

rhe results were obtained with a straight step at angles of 
sweepback from 0° to 57° on the 6-in.-wide sidewall of a blow 
down tunnel operating at JJ = 2.36 and at a Reynolds Number o 
9 & 10° per in 

Schlieren observation indicated a boundary-layer thickness of 
about 0.25 in. in the empty test section, and a step height of 0.31 
in. was chosen to ensure that, in the absence of sweepback, free 
separation would be obtained—that is, separation with pressure 
distribution independent of the agency provoking it With 
sweepback, the wall pressure distribution along a line normal to 
the separation line was characteristically similar in shape to that 
obtained at zero sweepback, and the separation line was just as 
far ahead of the step, so that it would appear reasonable to con 
clude that free separation was again achieved 

\t each sweepback angle the streak patterns obtained from 
applying the China film technique showed no variation along any 
Straight line parallel to the separation line for at least 1 in. on 
either side of the line of wall pressure transverse. It was therefore 
concluded that there were no changes in the flow in the direction 
parallel to the separation line 

Fig. | shows the measured ratio of peak static pressure to free 
stream static pressure plotted against the component normal to 
the separation line of the main-stream Mach Number upstream 
of separation. The Figure also includes a point from an experi 
ment where a state approaching freet separation with sweepback 


* This investigation was conducted under a research grant from the Uni 
versity of Sydney during leave of absence granted the author by Aus 
tralian National Airways rhe guidance of Dr. R. E. Meyer is gratefully 
acknowledged 

’ The China film pattern showed a separated region about six initial 


oOundar ] 
boundary-layer thicknesses wide 
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Fic. | Effect of normal Mach Number eparati peak 
pressure ratic 
was obtained on the sidewall by means of the nose shock from 


wedge on the tunnel floor; the peak pressure in this case was taken 
to be the pressure at the characteristic inflection point of the 
pressure distribution between the separation and reattachment 
processes 

Comparison with results of previous investigations for the 
two-dimensional free separation shows that, within the margin of 
uncertainty of those results, sweepback principle may be 
enunciated-—namely, that the peak-pressure ratio for turbulent, 
free separation depends only on the Mach Number component 
normal to the separation ling 

This conclusion has been shown to be plausible by extending to 
three dimensions the model used successfully by Crocco and 
Probstein® to estimate the peak-pressure ratio at turbulent separa 
tion in two-dimensional supersonic flow Phis will be published 


in a later paper 
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Some General Properties of the Dynamic 
Amplification Spectra 


Y. C. Fung 
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ws A MORE OR LESS rapidly applied force acts on an 
elastic structure, the concept of dynamic amplification 1s 
of great importance. In the following we shall consider simple 
pulses of the nature of the landing impact loads; they rise mono 
tonically from zero to a peak value at time /,,, then fall off back to 


zero at time ¢ A plot of the maximum response versus fre 
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quency, of a variable frequency single degree of freedom linear we obtain from Eqs. (4), (5), and (7 
oscillator to such a pulse, normalized in some specific manner, , j ‘ 
R(t = W) Mo Tt uw Ft 1/2 Mi-/ pM t+ (« Lu 
is called a ‘‘shock spectrum Such a shock spectrum can be ob 
. . . >) 3 > ? j ‘ t 
tained experimentally, for example, by multifrequency reed gages 1/2) ( 1(e/3 Ee") M <€( MiMe/ i EU Mi") Mom) 5 & 
Since the idea of shock spectrum has been used broadly in specify 
. F m ‘ ‘ 14) Mi bh ¢ T €|( pfo/p bh bh - 
ing the nature of shocks in many engineering design problems 
° . ) 2 » 
see references 1 and 2 where a survey of literature is given), it piler/e 1/6 Ma) Me Mike) Me 
seems of interest to study the general features of such a spectrum } /D)@2 1/2 1/3) (g ( 
that are not particularly sensitive to the specific details of the 
° - /9? 5 
pulse history = exp|—e 2)) iu » (mr? /2y 
Let 3(t) be a pulse for which peo /2 € 1/2 m 1/2) (r?/p ; 
: , {) j< {) , / Ml ‘ 
rf for < and f ; « /6 e*/2 € Siu €or «€ Mi ple / 
QO for tin (0,!¢ 
(/Ode l r)€ vi iv Sa 
mlaXx y(t ¥(t l 
0 Eq. (6) shows again tp, > 7/2wasw—>0O. Eq. (11 
: . : : , the slope of the spectrum at the origin (@ 0 ] 
Phen the positive dynamic amplification spectrum with a small ; IS SPCCUMEAY GE CBS OFIEIE Lt is equal t 
i i i impuls multiplied by the nur cal or w/t 
damping parameter, e€, 0< «qi is defined as Pe Ne tiplied ry the umerical factor : 
bd | 
—— lim |d/d(at,,/x ) 7 ¥(r) dr 
. 
j | >) J 0 
man w jrje * sin w(t r jdt Z Se 
0 | J0 { It is clear that the curvature of a mplificati p 
; ae ; the origin (« 0) always vanishes and that tl igh 
Phe negative amplification spectrum is obtained by the same ; 
: ae terms can be expressed in terms of the successive mo 
expressions through replacing max by min hese spectra are : : 
; : : : . pulse about the origin 
considered as functions of the dimensionless parameter w TT 
P } o tl } ‘ i é 2) one i 1 sul 
Expanding the sine term in Eq. (2), one can rewrite the result Ren INvVEesk Peopren 
is 
; Po the inverse problem of determining 3 from the 1 
woe 4 Sin (@ so ad the following question may be askec Ca 1 
where mined from I xp ding into power 
a ind « paring coefficients in Eq 1), one obtains, i 
R we (| (7 )e®? cos wr dr + € U 
\Jo 
z yu? u d/d ) 
¥( 7 ie sin dr) ( } “ Pa - 1/3 d | 
\/7 0 
. t In general, for G; 1.2, 
— 2 eo 7 
ircta K( re in wrdr 
. Mo; f ~ LL LU 
es 0 S z NS 
— -\1O rit y a) ) +\"9 ) 
. j 0 < ~/) é é I é “ 
| ¥( res’? cose rar 5) | 2 | ; 
. o 
Here all t derivatives of or must be eva 
PROPERTIES OF THE (AMPLIFICATION SPECTRUM AT SMALL VALUES t< 0. It is seen that any w equations of the type 
OF FREQUENCY volve re than 7 s; thus the possibility of determin 
Now it can be shown that when w/ is sufficiently small, the inverse becor mcertain The same is true for t 
i R , : 1 ' . 1 oe } a sol al 
tbsolute maximum of g occurs at t > ¢ i.e., after the pulse has U<« However, it is yet unknown ether detained 
fallen back to zero. For instance, for a trapezoidal pulse thi mation oO1 it larger values of w v uniquely defi 
occurs when wt, / 7 1/3. That this is in getieral the case can be Ver St Phe negative spectrum, , and the acer 
een from Eq $) and (5) which show that 6 —~ O when — () spectra, , ete., are much more sensitive to the pulse 
ind that. for small n in ting theoretical question is, therefore, the possi 
" dete ling the inverse from such spectra 
R—>w ¥( Tr dr 
( : = 
adie LARGE VALUES OF FREQUENCY 
which increases monotonically with increasing ¢ if F(¢) is a simple ' ' ' 
; ; : : , ” 2 pie In what follows we will consider the no-damping cas¢ 
pulse defined by Eq. (1 Hence, the peak of g(t) is controlled, . 
é INE 
as in Eq. (3), by sin wt when wis small. If ¢, denotes the time at 
. 
which the first peak of g(?f) 1s reached, we have t, —@ 7/(2w), : 
: ’ ‘ * ‘ y T SIll wr dT 
which can be made fy by choosing a sufficiently small w no 
Since ¥(7 Qfor rt to, R and @ become constant for ¢ > f¢ . 
Hence, when wl,,/m is sufficiently small so that the first maxi r dsx(t r)/d r)] « r 
mum of occurs at tf > fy, Eq. (3) shows that the peak time, ¢), vv 
ind the peak value of ire given by co lim t)/q =f 
wl o = 7/2 6 
‘ ) a whicl hows that, if the applied load v iries slowly in cor 
Exp } el(7/2) 4+ ; Rt r - ; : : 5 oe 
with the natural vibration frequeney of the oscillator, 
respectively, where @» and R(fy) are the values of @ and R(t sponse tends to be identical with the pulss 
evaluated at Expanding ¢€*7 cos wr and e€*7 sin wr into \symptotie formulas for large values of w can be d 
lavlor’s series in w and defining the successive “moments” of the i method given by Willi Consider 1 following 
pulse ¥(7) about the origin as u which is related to the last term in Eq. (14 
. , 
m r 3(7r) dr, k = 0, 1, 2, be! Ila f(t r)< - @ dr 
WV 
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I r) into a Tavlor’s series with remainder, 
. 
Da — yr 0) 4 
| —_— 
J0 0 
T n'| f Af coos wre sales dr 17 
whe f 1. Expanding the following integral into a Taylor’s 
. 
COS WT ¢ = dr = Vit,a IS 
/J0 
Wit,a)= > 1,(t)a 19 
P a 
) 0 
then A,(t) = (1/r!) (0° p/0a" )a 20 


By differentiating Eq. (18) 7 times with respect to a, we obtain 
t tr)’ coswre” = "1 dr = ( ] Oy 0a 21 
J 0 


Substituting into Eq. (17) and passing to the limit a = 0, we have 


£ 


[(0, w a A 0 l ae | l 6t) A yy 4 
pe 
? 0 
N W 
*/ 
yf a cosare © dr = 
27 0 
w/(w? + a sin wt + (a/w) cos wf a/a - 
He ice 
1, = (1/w) sin wf 
1, = (1/w?) (cos wt — 1 = 
1, = (1/w*) (wl sin wt 
i; = 1 /w*) [(w?t?/2!) + cos wt l 
et From Eqs. (22), (23), and (14), we obtain the asymptotic 
formula 
F(t 1/w) sin wt #(0) + (1/w?) X 
cos wl 1) +(0 t F(t 4+- (1 /w*) sin wt 3 61), 
(O<@e<1 4 
t)/qo = F(t (1/w) sin wt F(O) + 
1/w?) |(cos wt 1) +10 t ¥°)(0) (12/2!) F°2)(6t 4. 
(1/w*) sin wt (0) (1/w*) (cos wh 1) ¥‘*)(6t), 
(Of 6< 1 (25 
et rhe indeterminacy of the factor @ can be removed if F(¢ 
can be approximated by a polynomial of degree n. Then F 


t) vanishes, and an expansion of the following form may be 


obtained: 


= F(t (1/w) sin wt F(0) + 


n—2 
1 /w? cos wt 1) #(0) Yt r!) Fx” TAO + 
; 1 


(1/w*) sin wt F (0) 4+ 26 
Eq 26) shows that the time tp at which Qmar is reached ap- 


proaches tm when w is large: 


tpt, 2% 
Hence, Imar/do = 1 ¥(0)/w] sin why + O(1/w? 28 


rhis shows that the wavy curve of the spectrum that converges 


to 1 as w increases has a “period” 


wlm = 2 29 
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Simplified Equations for Calculating Local and 
Total Heat Flux to Nonisothermal Surfaces’ 
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SYMBOLS 


width of surface, ft 
local heat-transfer coefficient, g } Bt I rf | 


thermal conductivity, B.t.u. ‘hour ft. | 


exponent in equation describing velocity variation 
P Prandtl Number, P mr 
q local heat flow per unit area and unit time, Bt I ft 
@) total heat flow per unit time, B t.u_ hour 
Re local Reynolds Number based on R 
absolute temperature, R 
velocity in x direction 
coordinate measured along the surface indicating ation 
v difference between wall and recovery temperature 
coordinate measured along surface ndicatin cation of ste 


wall temperature, ft 


Sudscrip:s 


= isothermal wall 

QO =} value at leading edge 
refers to adiabatic wall or recovery temperature 
refers to conditions at the surface 


boundar 


refers to conditions at location x outside the 


INTRODUCTION 


U" TO THE PRESENT TIME, the major efforts in the study of 
heat transfer to external surfaces have been directed to 

constant wall-temperature condition. Unfortunately, in many 
practical cases, the wall temperature is not constant but is vary 
ing along the surface The question then arises of how much the 
heat transfer in this case deviates from that on a constant wall 
temperature surface. Contributions to the clarification of this 
problem have been made by Fage and Falkner, Rubesin and 
Chapman, Bond, Eckert, Seban, Tribus, and others. An ex 
cellent review of such solutions is available rhe results of 
these analyses are somewhat difficult to apply to a completely 
arbitrary wall-temperature distribution. Consequently, there 
has been a demand for interpretation of these results in such a 
form that may be readily utilized by design engineers rhe pri 
mary purpose of this note is to present such a simplified method 
for calculating the heat transfer to surfaces for either completely 
laminar or turbulent boundary-layer flows with and without 
pressure gradients over surfaces with arbitrarily prescribed non 


uniform temperature 


Basic RESULTS 


It has previously been shown! that the local heat rate g(x) for a 
constant property boundary layer with arbitrarily varying wall 


temperature is rigorously given by the following expression 


Px 
Gx) = | h(x, ¢) (dd,,/d¢) dé 4 
J0 
> 
df a) : l 
In this equation, 3, = (7 T.) indicates the difference between 
the wall temperature and the recovery temperature, both taken 


at position ¢; A(x, ¢) is the heat-transfer coeflicient existing at 
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Fic. 1 Laminar and turbulent heat transfer from surface with 
pressure gradient and variable wall temperature 


location x when the surface is unheated upstream of position ¢ 
and heated to a constant temperature downstream of position 
¢. The integral portion is evaluated without taking into con 
sideration any discontinuities in the wall temperature. Such dis 
continuities are accounted for by the summation expression, with 
each discontinuity giving a separate contribution to the local 
heat flow. The symbol 3,(¢;*) indicates the temperature dif- 
ference just upstream and #,,(¢;~) the value just downstream of 
the temperature discontinuity. 

The required expressions for the local heat-transfer coefficient 
h(x, ¢) for a step in the wall temperature as obtained by approxi 
mate methods are available.! They will be referred to the heat 


transfer coefficient f;,, of a surface with locally constant tem 


perature. For laminar flow, 
h(x, £)/hio = [1 — (¢/x)* 3 (2 
The value of @ is dependent on the pressure gradient For 
wedge-tyvpe flows described by u., = cx”, the following results 
apply: 
m O 1/4 1/2 3/4 1 


In the case of turbulent flow with constant velocity over a flat 
plate, two approximate expressions are available. Both will be 
reported here since experimental data are not sufficient to decide 


which one is better 


Seban: h(x, ¢)/h = [1 i<c/2Fr 1/9 3 
Rubesin: h(x, ¢)/h = [] (¢/x )39/40] —7/39 | 


When the only temperature discontinuity occurs at the lead- 
ing edge, the expression for the local heat rate is 


eX - 
hi x» ¢) dd re 
\ =i vy, + : dt (2 
vl 


i d¢ 

A rather significant result is readily obtained if one assumes a 
wall-temperature distribution given by the expression J,, = Ax 
For this case the integral in Eq. (5) may be evaluated rhe re 
sults are shown on Fig. 1. It is seen that the influence of the 
nonisothermal wall conditions is more pronounced in laminar than 
in turbulent flow It may, however, be significant even in 
turbulent flow if severe temperature gradients are present. <A 
pressure gradient slightly decreases the influence of the variable 
wall temperature as demonstrated in Fig. 1 for laminar wedge 
flow 

The present note utilizes a method suggested in reference 2 
The actual wall-temperature distribution is approximated by a 
series of sloping segments as shown in Fig. 2, thereby making 
dJ,,./dé a constant value over any interval Ax and allowing its 


removal from the integral sign. Eq. (5) then reads 


TICAL SCIENCES—JULY, 195 





——_+ 


| | 
UF Tw > Tr | | 








| | | 
| | | : 
| AX» 
) | | , 
| | g 
1 1 1 l 
Xn-1 Xn 
Xx ( 


Fic. 2. Laminar sloping segment approximation for temperatur 
distribution | 


general form of the integrals in this equation may be 


The 


written as follows: H 


F h(x, ¢ Pane Bis of ¢ 
dj = x | d 
. h J0 h \ 


; Stré 
Since the expressions for A(x, ¢)/h;.. are known functions of of 
¢/x) for laminar and turbulent flow, it is possible to evaluat 

the integrals once and for all 
The values of the integral on the right side have been obtaine 

ind tabulated in reference 3 for values of (x;/x) from O to Lt 

steps of 0.01. A further simplification results from the fact tl 

the tabulated values of the integral can be approximated wit 


good accuracy by the following expression: 


5 


+ Me, 2 c \ \ 
d =A + B S 
J 0 h "4 v \ 


j t 


When this relationship is substituted into Eq. (6) under t 
assumption that the length x is subdivided into » equal intervals 
i.ec., Ax = constant—we arrive at the final relatively simple ex 


pression for the local heat flux g(a 


a(x = h ,od + A(d v + B( Ax /a x 
) I ‘ 


v 2d 29 28 + 2n 1)’ { 


[he dimensionless parameters in this equation and in the equ 
tion hiss = C(k/x)Rer" Pi for the isothermal heat-transfef 
coefficient are collected in Table 1 

It is estimated that Eq. (9) will vield results which are within 


5 per cent of those obtained by using the more exact value of the 
integral as reported in reference 3 
TABLE 1 I 
m A B ( 
Laminar 0 0.895 0.690 0.332 1/2 
1/4 0.837 0.635 0.412 1/2 
1/2 0.840 0.572 0.469 1/2 
] 0.792 0.538 0. 560 1/2 
rurbulent 
Seban 0 0.99] 0.117 0.0296 0) 
Rubesin 0 0.967 0.256 0.0296 Ss 
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~ TABLE 2 


Lat () 0). 696 0.432 0.664 ifs ® 


Tur nt : / ai 

po Sebat 0. 982 0.478 0.0370 0.8 | * 1 Yad ‘ae 
: R ; 1.018 0.475 0.0370 OLS C3 t 

| | © CNT OAL 












. 4° rmx a ee eS 
| ; ° 1 
UATION OF THE TotraL HEAT FLux Q(a r-§- ‘ 
- 22 ~ 
It 1 be desirable to evaluate the total heat transfer up toa - "s - 
giv ition along the surface This may be accomplished Fic. 1 Sketch of specimen 
inte ting the local heat flux as determined above A direct 
| pproa in be made for the flat plate by performing this inte 
| gration with the same basic procedur¢ is outlined in the preceding 
— secti Phe resulting expression for the total heat flow up to a 16 
sition xX, mav be expressed S330 IN/MIN 
° 
A) ] vb} (od 4 1'(9 J B'( Ax/x x 4 T T 
) 2 DY L (9x ' (10 
mperatur 1 ot Uw, 2n l v { l 
12} 


iensionless parameters in this equation and in the expres 





nh C'(k/x)Re"Pr are contained in Table 2 


| of 6,600 IN/MIN ||. 


results obtained with Eq. (10) are in substantial agree- 





ULTIMATE NET STRESS —PSix 1073 














1 1 : . 4 co 
nent (within 1 per cent) with those obtained with the tabulated peo cea 
1 e 
integrals for both laminar and turbulent flow 8 +++++ *, #8 D2 ~/min 
o 
hd 
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Specimens were prepared by drying for 48 


Stress Concentration in Plexiglas as a Function the response needed 
cuions of of Loading Rate hours at 122°F. followed by at least 5 hours in a desiccator at 
evaluat room temperature. Room temperature was maintained between 


74. and 78°F. for all tests 

















Klir 
obtain Assistant Professor, Division of Engineering Sciences, Purdue \ total of 70 specimens designed as shown in Fig. 1 were used by 
Oto li ¢ y. Lafayette, Ind Ely of Redstone Arsenal to establish the upper curve of Fig. 2 
fact that April 4, 1957 For the present study, specimens were modified from the design 
ted wit shown in Fig. | by drilling a centrally located hole parallel to the 
SUMMARY 0.25-in. dimension. Specimens with hole diameters of ‘ 
; ind in. were studied. The net rupture stresses ‘ for 
dy made of stre concentration effect is functions of loading 

| . , | hol nl fon lo yortion 
S rat Pols methy! methacrylate is chosen as the working medium The the specimens with in. holes are plotted on the lower portion 
nditions used are uniform crosshead velocities ranging from 0.02 of Fig. 2 as functions of the rupture time rhe net rupture stres 
nder tl saniertbsclae  phiara sis Aes etcsaeliiand snes It is shown that failure to is considered to be the quotient of the rupture load and the net 

. le ence of the stress concentration on the loading rate may 1 1 1 
l area measured at the cross section through the rie Phe lower 

ervals de 
mple ex curve of Fig. 2 represents a reasonable fit for the experimental 
INTRODUCTION data 
r of the loading rate on the ultimate strength an Discuss 
' (9 dulus of elasticity of plastics has been the subject of con 

are os ; sing curves similar to Fig. 2. a strength-reduction factor K 
iderable research. The results of this research cannot be in Using curves similar to Fig. 2 ength-reaucth ctor 4 
ne cq can be calculated for the geometries included in this study. This 


igently used in designing plastic-like items such as rocket 


transfer = ” 


° - i Ss ] 7 >) 
r propellant grains until the dependence of stress concentra factor is defined to be 








loading rate is also established. This work supplies in K yr /S 

I t how stress concentration effects are influenced by , } | 
“aeeer where S, is the rupture stress for the material as determined from 
oe specimens with no stress concentration. This factor is to be dis 
EXPERIMENTAL WORK tinguished from the theoretical stress concentation factor Ay, a 

function of geometry only, which is defined b 
‘ ethyl methacrylate (Plexiglas IA, UVA) was tested in : , 
, vitl even different crosshead velocities ranging from A > ) ~ 
1/2 <1. per min. to ipproximately 6,600 in per min. An Instron where S.. is the maximum stress in the stressed body and 5S, is 
1/2 s used to obtain velocities from 0.02 to 2 in. per min the average stress over the net cross section of the stressed body 
1/2 rapid rate machine designed at Redstone Arsenal’ was used for \lthough the theoretical stress concentration factor 1s of con 
1/2 thre er velocities. Although this machine is not siderable value, the designer is usually more interested in the 
t crosshead velocity machine, the crosshead velocities strength-reduction factor A For a given combination of hol« 
2 : rmetr Research Engineer, Research Laboratori Orda=n M size and loading rate, the value of K,. can be obtained with the 


Redstone Arsenal iid of the curves similar to Fig. 2. For example, to find A, for 
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TABLE 1 
Stress Concentration Values 


the 5-in) per min. crosshead velocity and the !/s-in. diameter 
hole, the value of 11,700 psi for S, is read from the top curve of 
Fig. 2 while the value of 7,950 psi for S,, is read from the lower 
curve. Thus 
K, = 11,700/7,950 = 1.47 3 
Table 1 includes the calculated values of A, for each combination 
of stress concentration and loading rate 
Values of the theoretical stress concentration factor? A, for the 
hole sizes studied are If, for 


geometry, A, is assumed to be constant for all values of maximum 


included in Table 1. a particular 


normal stress S,,, then it is permissible to write 


Bie 2 S| Ses 


where S,,,; is the maximum normal stress at rupture. If the fur 
thur assumption is made that the maximum normal stress theory 
of failure accurately predicts the behavior of the material, then 


S,, and the theoretical stress concentration factor becomes 


Snr = 


Ky, = S,/Sue a 


which is the definition of A Since the experimentally obtained 
values of A, do not e- ual A 
in arriving at Eq (5) are not accurate 
tions (that K 
tic behavior 
perience sufficiently high strain rates in the vicinity of the concen 


it follows that the assumptions made 
The first of these assump 
is constant) is accurate if the material exhibits elas 


All of the specimens with stress concentrations ex 


tration so that the assumption of a linear stress-strain relation is 
reasonable. It is felt, therefore, that the first assumption used in 


arriving at Eq. (5) isa good assumption. Thus, the second of the 
above assumptions (concerning the theory of failure) is primarily 
responsible for the difference A, — A,. This difference can 
therefore be considered as a measure of the inability of the 
maximum normal stress theory of failure to predict the behavior 
of the material. The values of A, — A, 


As the loading rate is increased, the material loses ductility 


are included in Table 1 


Since the maximum normal stress theory of failure is most ac 
curate for brittle materials, it might be expected that A, — A, 
would decrease as the loading rate increases. This is the case as 
can be seen from the tabulated values of A, — K, 


CONCLUSIONS 


Failure to recognize the dependence of A, on the loading rate 
may result in serious underdesign. Consider the case where the 
tress distiibution around a point in a propellant grain can be con- 
sidered comparable to the distribution around a !/s-in. hole in a 
tensile specimen loaded with a crosshead velocity of 330 in 
Using the value of A, for a commonly employed crosshead 


per 
min 


rate, 0.02 in. per min., gives a design rupture stress of 


S = 14,900 12,100 psi 6) 


1.23 = 


Using the value of A, which should be associated with a crosshead 
velocity of 330 in. per min. gives a design rupture stress of 
S= 167 = 


14,900 8,920 psi 7) 


For this example, failure to recognize the dependence of A, on the 
loading rate would cause one to predict a design stress 36 per cent 
higher than the correct value. 


rrcCAL SEILENCES 


The results of this study show that the strer 


AK, increases as much as 58 per cent over the range 
considered 
The factor KA; — A, deere 


This is an indication that the maximum normal stress t 


ises with increasing 


failure becomes more accurate as the loading rate is incré¢ 
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Errata—‘‘Vibration of Clamped 
Parallelogrammic Isotropic Flat Plates’’* 


M. Hasegawa 
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4 ipo: TERM + Ow?*/Oy' should be added 
a(Otw/OxOv*) within the bracket in Eq. (1 


p should read p in Eq. (38 


after 


p? should read p* in Eq. (6) and in the first line of Eqs 
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r | SHE FOLLOWING two errors exist in Table 1: 


(1) The value of w*Zg for w = .40 should be 


— .68350 + .31076 i 


(2) The value of w?.\/g for w = .20 should be 


60826 049698 i 
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